An accelerated stochastic vortex structure method for particle collision and agglomeration in homogeneous turbulence

Modeling the response of interacting particles, droplets, or bubbles to subgrid-scale fluctuations in turbulent flows is a long-standing challenge in multiphase flow simulations using the Reynolds-Averaged Navier-Stokes approach. The problem also arises for large-eddy simulation for sufficiently small values of the Kolmogorov-scale particle Stokes number. This paper expands on a recently proposed stochastic vortex structure (SVS) method for modeling of turbulence fluctuations for colliding or otherwise interacting particles. An accelerated version of the SVS method was developed using the fast multipole expansion and local Taylor expansion approach, which reduces computation speed by two orders of magnitude compared to the original SVS method. Detailed comparisons are presented showing close agreement of the energy spectrum and probability density functions of various fields between the SVS computational model, direct numerical simulation (DNS) results, and various theoretical and experimental results fou...


I. INTRODUCTION
Particle collision and agglomeration play an important role in a wide range of turbulent flows applications involving small particles or droplets. Droplet collision is a key element to cloud formation and precipitation development (Devenish et al., 2012). Particle agglomeration is particularly important in aerosol flow problems, such as fly ash collection from combustion processes (Xu et al., 2010), flame-synthesis of nanoparticles (Zhang et al., 2012), cyclone particle separators (Paiva et al., 2010), and snow crystal formation (Kajikawa et al., 2000), for which adhesive particles have Stokes numbers sufficiently close to unity that they display significant drift relative to the fluid trajectories. Agglomerate formation is preceded by particle collision, where the particle collision rate is controlled either by the fluctuating turbulent shear flow (for smaller size particles) or by particle inertia (for larger particles). The fluctuating turbulent shear stress also controls agglomerate breakup (Serra et al., 1997 andHigashitani et al., 2001). Over long time, the distribution of particle agglomerate sizes is determined by a balance between influences increasing collision rate and influences enhancing agglomerate breakup.
A wide variety of turbulence models have been developed using the Reynolds-averaged Navier-Stokes (RANS) approach, ranging from the popular two-equation models, such as k − ε and k − ω, to full Reynolds stress models. RANS models yield numerical predictions for the mean turbulent a) Author to whom correspondence should be addressed. Electronic mail: jmarsha1@uvm.edu. Telephone: 1 (802) 656-3826. velocity field as well as for certain averaged quantities associated with the Reynolds stress tensor. However, additional modeling is required for RANS simulations to account for the role of turbulent fluctuations on particle transport. A similar need for subgrid-scale modeling of turbulent fluctuations arises for large eddy simulations (LES) when the Kolmogorov-scale Stokes number is less than a critical value of about three (Jin et al., 2010).
While numerous effective methods are available to simulate the effect of subgrid-scale fluctuations for transport of non-interacting particles (e.g., Wilson and Sawford, 1996;Loth, 2000;Minier et al., 2014;and Pope, 2011), turbulent subgrid-scale simulation for interacting particles remains an unresolved modeling challenge. There are a number of reasons why subgrid-scale modeling for interacting particles poses difficulties. First, the mechanics of interacting particles depends sensitively on the distance between the particles at small values of separation. Particles that are sufficiently close to each other experience highly correlated fluid velocities induced by the nearby turbulent eddies. Models which employ independent (uncorrelated) stochastic forcing at each particle consequently cannot be used for interacting particles. Second, particle collision and adhesion processes occur on very small time scales, which makes the numerical simulation of colliding and adhesive particles numerically stiff. This is particularly a problem for simulations using the soft-sphere discrete element method (DEM), which is usually necessary for dealing with particle agglomerates that form upon collision of adhesive particles. Consequently, small time steps must be taken for the particle transport and the subgrid-scale turbulent fluctuation modeling must be sufficiently fast for the computation to be manageable. Third, the eddy structures of the turbulent flow play an important role both in dispersing particles and in inducing clustering in the region in-between the eddies (Squires and Eaton, 1991;Bec et al., 2007;Grits et al., 2006;and Falkovich and Pumir, 2004). Eddy-induced particle clustering leads to formation of regions of high particle concentration within the turbulence, which dramatically increases particle collision rate and agglomerate sizes (Sundaram and Collins, 1997;Zaichik et al., 2006;and Reade and Collins, 2000). Particle preferential concentration has particularly interesting consequences in bidisperse flows involving particles that are both heavier and lighter than the fluid, such as heavy particles and bubbles in a liquid (Fayed, 2013 andFayed andRagab, 2013), for which case the heavy particles cluster in the high shear regions in-between the eddies and the light particles (bubbles) cluster within the turbulent eddies. As a consequence of the issues of computation time and preferential concentration, many of the synthetic turbulence approaches that have been developed for reconstruction of initial or inlet conditions in large-eddy simulations (Kraichnan, 1970;Smirnov et al., 2001;Tabor and Baba-Ahmadi, 2010;and Lund et al., 1998) are not useful for subgrid-scale modeling of flows with interacting particles.
Clustering of non-adhesive particles in turbulent flows is largely due to inertial particle drift across curved fluid streamlines associated with the presence of turbulent eddies (Squires and Eaton, 1991). A vortex structure representation of the turbulent flow consequently presents a natural approach for capturing this effect in the turbulence model. Of course, vortex-based structural models have long been discussed in the turbulence flow literature. Notable among these are Townsend's (1951) model of homogeneous turbulence as a collection of Burger's vortices and Lundgren's (1982) spiral vortex model of turbulence. The scaling and structure of coherent vortices were examined by Jiménez et al. (1993) in homogeneous turbulence based on results of high-resolution direct numerical simulations (DNSs) and by Berlin et al. (1996) in a turbulent shear flow using experiments with low-temperature helium gas. Both studies found that the vorticity field for low Reynolds number turbulence is dominated by a set of strong, coherent vortex structures of finite length and with tubular shape, surrounded by a sea of weak random (non-coherent) vorticities. The length and core radius of the coherent vortices were found to scale with the integral length scale and the Kolmogorov length scale, respectively, and the vortex strength was found to scale with the square root of the microscale Reynolds number. Analysis of the Townsend and Lundgren models was given by Pullin and Saffman (1993) and Saffman (1997), who derive an expression for the energy spectrum and other measures for isotropic turbulence. Kivotides and Leonard (2003) report results of a computation in which homogeneous turbulence is represented by a set of finite-length vortex structures, and show that this system generates an energy spectrum that satisfies the Kolmogorov k −5/3 scaling in the turbulence inertial range. The effect of vortex straining on the energy spectrum of a group of randomly advected vortices is discussed by Malik and Vassilicos (1996). Hatakeyama and Kambe (1997) demonstrate good agreement for structure functions of homogeneous turbulence between those generated by a group of random strained Burgers vortices and the classical Kolmogorov theory. Use of vortex models to generate accurate PDF curves for velocity increment, acceleration, and vorticity is discussed by Min et al. (1996), Wilczek et al. (2008), and Wilczek and Friedrich (2009).
A first step toward the use of a vortex structure model for turbulent particle transport was made by Ayyalasomayajula et al. (2008), who proposed a model in which the turbulent eddies are represented by a two-dimensional vortex array and a stochastic algorithm is used to vary the strength of each vortex in time. Although extremely simple, this model was shown to yield reasonable results for particle acceleration statistics and clustering. A three-dimensional stochastic vortex structure (SVS) model was proposed by Sala and Marshall (2013), in which the turbulent vorticity field is approximated by a set of finite-length, fixed vortex structures which are randomly positioned and oriented in the flow field. Predictions of the SVS model for turbulence energy spectrum and particle collision rate were found to be in close agreement with DNS predictions. However, the original SVS method was rather slow and only considered transport and collision of non-adhesive particles.
The current paper extends the SVS model proposed by Sala and Marshall (2013) in two respects: (1) a variation of the fast multipole method FMM and local Taylor expansions are used to dramatically accelerate the SVS computations and (2) the performance of the SVS method is examined for prediction of turbulent agglomeration of adhesive particles. Successful simulation of turbulent agglomeration requires both that the particle collision model is accurately simulated by SVS and also that the fluctuating turbulent shear stress responsible for agglomerate breakup and erosion is accurately predicted. We also report more extensive comparisons with DNS data, as well as detailed sensitivity testing of the SVS model results to various input parameters. The basic SVS model is described in Section II. In Section III, a fast multipole method is developed for computing the velocity field induced by the vortex structures, which is found to yield nearly two orders of magnitude increase in computational speed compared to direct velocity computation. Sections IV-VI present different types of validation and sensitivity tests for the SVS model. Section IV examines measures of the turbulent flow field. Section V examines prediction of collision rate for non-adhesive particles and Section VI examines use of SVS for prediction of turbulent agglomeration with adhesive particles. Conclusions are given in Section VII.

II. STOCHASTIC VORTEX STRUCTURE METHOD
Particle collisions in turbulent flows depend primarily on the eddy Stokes number, which can be written as a function of eddy size ℓ as where d and m are the particle diameter and mass, respectively, and µ is the fluid viscosity. In the inertial range, the characteristic velocity u ℓ of eddies of size ℓ varies with turbulence dissipation rate per unit mass ε as u ℓ ∼ (εℓ) 1/3 (Frisch, 1995). Since the dissipation rate is approximately independent of scale in the inertial range, the Stokes number varies with ℓ approximately as St ℓ ∼ ℓ −2/3 . Particles are largely transported by the fluid flow for eddies where St ℓ ≪ 1 and the particle inertia filters out the turbulence fluctuations for eddies where St ℓ ≫ 1 (Ayyalasomayajula et al., 2008 andMarshall andLi, 2014). In-between these extremes, there exists an eddy size ℓ for which St ℓ = O(1), in which the particles are thrown out of the turbulent eddies and collect in high-concentration sheets in the interstitial region between the eddies. The stochastic vortex structure (SVS) model approximates the turbulent vorticity field by a collection of vortex structures placed and oriented randomly in the flow field. In the simplest version of the SVS method, the vortex structures all have the same finite length L, core radius δ, and strength Γ, although a multiscale version of the SVS model has also been developed. The vortex length L is assumed in the current paper to be of the order of magnitude of the turbulence integral length scale ℓ 0 = 0.5 u 0 3 /ε, where u 0 is the turbulence root-mean-square velocity. Based on the well-established observation that strain rate in the inertial range scales as u 0 /λ (Frisch, 1995), Kambe and Hatakeyama (2000) used a scaling analysis to derive an approximation for vortex core radius as δ 3.9η, where λ = u 0 (15ν/ε) 1/2 is the Taylor microscale, η = (ν 3 /ε) 1/4 is the Kolmogorov length scale, and ν is the kinematic viscosity. This estimate is in good agreement with experimental and numerical results (Jiménez et al., 1993 andBerlin et al., 1996). The current paper uses a somewhat larger assumption δ = 8η for vortex core radius in order to ensure sufficient number of grid points to adequately resolve the velocity gradient across the vortex cores. Each vortex structure has a lifetime T V which is assumed to be proportional to the integral time scale, T ℓ = q/3ε, where q = 1.5 u 2 0 is the turbulent kinetic energy per unit mass. While the coherent vortices in a turbulent flow may in practice last significantly longer than T ℓ , the results of the model are not sensitive to value of T V . The initial age of the nth vortex structure, τ 0n , is specified as a random variable, where the ratio τ 0n /T V has a uniform distribution between 0 and 1. If t 0n denotes the time at which the vortex structure is initiated, then the current age of the vortex structure τ n (t) is given by When τ n (t) exceeds the specified lifespan T V , the vortex structure is removed and a new vortex structure is introduced with random position and orientation in the flow. The vortex structures induce a velocity field u which is computed using the accelerated method described in Section III. Each of the N V vortex structures is advected in time by moving the two endpoints of the vortex structure by solving where the index n identifies the vortex structure and i (=1, 2) identifies the endpoint of the structure under consideration. After moving the end points, the vortex length is reset to L. The centroid position x n and unit tangent vector λ n for each structure are then recomputed from the positions of the new endpoint locations.

III. ACCELERATED METHOD FOR VELOCITY CALCULATION
The stochastic vortex structures constitute a kinematic representation of the turbulent flow, which is intended to generate a synthetic fluctuating velocity field that exhibits similar statistical properties to the actual turbulent flow. The dynamics of the turbulent flow is simulated by whatever RANS model is used to compute parameters such as turbulent kinetic energy and dissipation rate, and not via the SVS model. With this point in mind, it is recalled that a divergence-free vorticity field ω can be generated from the vorticity ω * associated with a set of finite-length vortex tubes as where Substituting (4) into the Biot-Savart equation, where s ≡ |s| ≡ |x − x ′ |, and using Green's theorem, one can readily show that the ∇ζ term in (4) makes no contribution to the induced velocity field (see the Appendix).
For computation of particle transport, it is more efficient to compute the fluid velocity on a Cartesian grid covering the computational domain and then interpolate the velocity from the grid nodes onto the Lagrangian particles that move through the grid. This is particularly true when using a multiple time scale algorithm for particle transfer (Marshall, 2009), in which the fluid velocity is computed on a larger time step than that used for transport of either free or colliding particles. The computations in the current paper are performed using a 128 3 Cartesian grid to cover a cubic, triply periodic domain with side length 2π.
To accelerate the velocity computation, we utilize the combination of an optimized fast multipole method (FMM) for computing the velocity field induced by sufficiently distant vortex structures and a local Taylor expansion to reduce the number of points at which the Biot-Savart integral must be solved. The accelerated method is based on a partitioning of the computational domain into a tree family of boxes consisting of some number M levels, each of which covers all grid points in the domain. The first level (m = 1) consists of the entire grid and has only one box. The second level (m = 2) consists of 8 boxes, which are obtained by dividing the side length of each box in level 1 by a factor of two, as illustrated in Figure 1. This division process is repeated for subsequent levels, with the number of boxes in each level m increasing as 8 m−1 . The boxes associated with the highest level are called the small boxes of the box family.
The velocity is evaluated at each point of the Cartesian grid by solving for the contribution to the Biot-Savart integral (6) from all vortex structures in the computational domain, as well as from neighboring domains necessary to enforce the periodic boundary condition. In order to perform the computation efficiently, we first associate with each grid point a specific smallest box of the tree family in which the grid point is contained, which is called the target box of the grid point. The velocity within each target box is determined by integrating the Biot-Savart integral over the vortex structures contained within some set of boxes (called source boxes) that can be at any level of the box tree family, but where the set of source boxes is required to cover each vortex structure within the computational domain exactly once (i.e., the source boxes cannot overlap). Each target box interacts with each source box either directly or indirectly. In a direct interaction, the velocity induced by each vortex structure in the source box is evaluated individually on each grid point within the target box. In an indirect interaction, the induced velocity from all vortex structures within the source box is computed at the center of the target box at one time using a multipole expansion, and then the induced velocity is extrapolated onto the grid points within the target box by a local Taylor series expansion. Lists are compiled for each target box of source boxes with which the target box interacts directly and indirectly. The selection of source boxes and the box interaction lists were constructed using the optimized approach proposed by Marshall et al. (2000), which is based on an analytical error estimate for the multipole expansion derived by Salmon and Warren (1994).
FIG. 1. Image representing two levels of the box family used to cover the computational grid. The first level consists of the entire grid and the second level consists of the eight individual boxes numbered 1-8 in the image. An example is shown where box 1 is a source box (blue online) and box 7 is a target box (red online), where the vector pointing from the centroid of box 1 to that of box 7 is indicated by an arrow and denoted by r. The individual vortex structures contained within box 1 are represented by short line segments within the box.

A. Direct velocity computation-Interpolation from the data plane
For a source box that interacts directly with a given target box, the velocity induced by each vortex structure in the source box is computed at each grid point in the target box. The velocity computation is done by first pre-computing the velocity induced by a vortex structure of unit strength on the data plane, which is defined as the positive r-z plane relative to the axis of the vortex structure ( Figure 2). This computation is performed once at the beginning of the computation and the results are stored.
The induced velocity on the data plane is determined by computing the induced velocity normal to the r-z plane of a coordinate system that is local to a vortex structure of unit strength, where the vortex center is located at the origin of the local coordinate system. The velocity at each point of the grid used to cover the data plane is determined using a Gaussian vortex blob method (Marshall and Grant, 1996), where the number of vortex blobs, N b , used to discretize the vortex structure is set equal to N b = int( βL/δ), and where the Gaussian radius of the blob is set equal to the vortex structure radius δ and β is a blob overlap coefficient. If the centroid of the ith vortex blob is denoted by b i , i = 1, . . . , N b , the associated vorticity field is given by Here, the blob amplitude Ω i is given by and λ b is a unit vector tangent to the vortex structure axis. Substituting (7) into the Biot-Savart integral (6) yields the velocity field induced by the ith vortex blob at a target point x as where P(a, z) is the incomplete gamma function with limits P(a, 0) = 0 and P(a, ∞) = 1. When a = 3/2 and z = x 2 for some real variable x, a convenient expression for the incomplete gamma function in terms of the error function erf(x) can be written as (Abramowitz and Stegun, 1965) The velocity at any point x on the data plane is obtained by summing the velocity induced by all N b vortex blobs, as given by (9). At subsequent times, the induced velocity from a vortex structure m at grid point x is obtained by interpolation from the data plane. This interpolation is performed by centering the data plane at FIG. 2. Schematic diagram showing the interpolation procedure used for direct computation of the velocity induced by a vortex structure on a node of the grid cell. Here L is the length of the vortex structure and P identifies the inclined plane from which the induced azimuthal velocity v induced by the vortex is interpolated. the vortex structure centroid x m , and orienting the plane so that it passes through the target point x and is tangent to the vortex axis unit vector λ m , as illustrated in Figure 2. The grid cell in which lies the point x is obtained in the data plane by integer division and the velocity induced by the vortex structure is interpolated onto the target point and reoriented to lie in the global coordinate system, yielding a velocity contribution u ′ m on point x from vortex structure m. The periodic boundary condition is enforced by including velocity induced by vortex structures in one period on each side of the computational domain, resulting in 27N V total vortex structures if the entire computation is performed directly. The total direct velocity at a point x from the N dir vortex structures for all source boxes on the direct interaction list (including vortex structures in the side domains used to enforce periodic boundary conditions) is then given by Since the sum (11) must be computed for every grid point within the Cartesian grid, it is very time consuming if the summation is performed over all vortex structures in the computational domain and the neighboring periods of the computational domain. For this reason, the direct interaction list is restricted to only a small number of source boxes with centroids located sufficiently close to the centroid of the target box.

B. Indirect velocity computation-Multipole expansion
For a source box ℓ that interacts indirectly with the target box, the contribution of all vortex structures in box ℓ is evaluated at any point x in the target box using the multipole expansion (Greengard and Rokhlin, 1987) where r = x − ξ ℓ is the vector from the centroid ξ ℓ = ξ 1,ℓ e x + ξ 2,ℓ e y + ξ 3,ℓ e z of box ℓ to the point x.
The box moment I ℓ, mnk of box ℓ is defined by The box moments are evaluated by first computing the moment J mnk of a single vortex structure of unit strength aligned in thex-direction about the vortex centroid in a local coordinate system (x,ŷ,ẑ), which is given by whereω(x)λ is the vorticity field associated with the vortex structure and λ is a unit vector along the vortex axis. For a vortex structure of length L and core radius δ, we find Since the values of J mnk are isotropic (the same for all directions), it is not necessary to translate between the local coordinate system used to compute (14) and the global Cartesian coordinate system. The moment I ℓ, mnk of a box ℓ is obtained by summing over the moments J mnk of all of the N ℓ vortex structures in box ℓ, which have vortex strengths Γ i and centroid locations c i , giving Once the moments of all of the smallest size boxes are obtained using (16), the moments of higher-generation boxes are obtained from the translation formula where i denotes one of the eight offspring boxes of parent box ℓ.

C. Indirect velocity computation -Local Taylor series expansion
The derivative term in (12) depends on the location of the target point x. Since we compute the velocity at each point of a N 3 Cartesian grid, there are typically a large number of target points within a given box. The local expansion method accelerates the process of computing the indirect component of the velocity field by evaluating the velocity induced by a source box ℓ with centroid ξ ℓ at the centroidξ b of the target box b (defined as the smallest box containing the target point x), and then determining the velocity at each individual grid point x using a local Taylor series expansion of K(x − ξ ℓ ) = r/r 3 about the target box centerξ b , given by Substituting (18) into (12) and truncating the summation after P terms gives the contribution of source box ℓ to the velocity at grid point x as where

D. Example computations
A series of example computations were performed with 512 vortex structures on a 128 3 grid with B different levels of box division. The order p of terms in the multipole and local expansions was allowed to vary from p = 0 to a maximum of p = 2 for all of the remaining computations. The order of the interaction is set for each source-target box combination as a function of the distance d between the box centers. The critical separation distance for each order is specified as a function of the box size b at the highest level B, such that we use order  Table I for a case with critical distance coefficients α 0 = 4, α 1 = 3, and α 2 = 2. The table lists the computed value of turbulent kinetic energy (TKE) (a measure of accuracy), the CPU time, and the percentage of the total possible boxes placed on the direct list (averaged over all target points). The CPU time results are for single-processor calculations for ease of comparison. It is noted that some source boxes do not enclose any vortex structures, in which case the box is ignored and not placed on either the direct or indirect list. At the top of the table are data for a computation in which the velocity is computed TABLE I. Comparison of CPU time, percentage of the computation performed directly (in terms of number of boxes of the smallest size), and flow measures such as turbulent kinetic energy and enstrophy for the direct computation (set in bold type) and for indirect computations with four different levels of the boxing scheme used for the velocity acceleration method. The computations were performed with N v = 512 vortex structures, with critical distance coefficients α 0 = 4, α 1 = 3, and α 2 = 2. using only the direct interaction. For B ≤ 4, the TKE error for computations using the accelerated method is less than 1.5% of the direct computation, while the CPU time is reduced to less than 3% of that for the direct computation. The CPU time is reduced further for the case with B = 5 to about 1.5% of the direct computation time, but at the same time the TKE error increases to about 16%. The reason for this sudden increase in TKE error is that the box size b grows progressively smaller as B is increased, so that an increasingly large percentage of the computation is performed using the indirect approach. As discussed by Salmon and Warren (1994), the multipole expansion error increases in a nonlinear manner as the critical distance decreases. Based on the results in Table I, we selected to perform the remainder of the computations in the paper with B = 4 and α 0 , α 1 , α 2 = 4, 3, 2.

IV. ANALYSIS OF THE SVS SYNTHETIC TURBULENCE FIELD
The key parameters associated with the SVS method are the number of vortex structures N V in the computational domain, the strength of each vortex structure Γ, the vortex length L, and the vortex core radius δ. These parameters can be related to various measures of the turbulent flow field, such as the turbulent kinetic energy per unit mass E, the dissipation rate per unit mass ε, and the enstrophy per unit volume Ω, defined by where D i j are the components of the rate of deformation tensor, u and ω are the velocity and vorticity vectors, respectively, and V is the computational domain volume. For homogeneous, isotropic turbulence, the dissipation rate and the enstrophy are related by ε = 2νΩ. The enstrophy can be estimated using the expression for a Burgers vortex (Burgers, 1948) in a field with axial stretching rate c, in which the vorticity field has the form of a Gaussian, and the Gaussian radius is δ = 2 √ ν/c. For a system of N V Burgers vortices of length ℓ and strength Γ, the enstrophy is given by A theoretical expression for the energy spectral density e(k) in a system of N V Burgers vortices of length ℓ and strength Γ is given by Saffman (1997) as where k is the wavenumber magnitude. This expression is derived based on the assumption that the vortices do not interact with each other, so that the energy induced by each vortex can be added together to obtain the total system energy. Integrating over the wavenumber interval (k min , k max ) gives the turbulent kinetic energy as where E 1 (·) is the exponential integral function. A series of computations was performed in which the number of vortex structures in the computational domain was varied from 32 to 512, and the product N V Γ 2 varies from 0 to 4000. The velocity field is computed using the accelerated method described in Section III. The mean computed values of enstrophy Ω and turbulent kinetic energy E obtained from the definitions (22) are plotted as a function of N V Γ 2 in Figure 3(a). In both cases, the computational results collapse onto a single line, as predicted by (24) and (26). Since both enstrophy (and hence dissipation rate) and turbulent kinetic energy are proportional to the combination N V Γ 2 ℓ/V , the modeler is free to select N V based on an alternative criterion and then to set Γ to obtain the desired turbulent kinetic energy.
There is a slight variation in the computational values of turbulent kinetic energy and enstrophy depending on the randomly selected positions and orientations of the vortex structures. In order to characterize the amount of variation caused by the random character of the SVS algorithm, the turbulent kinetic energy and enstrophy calculations were repeated 10 times and the root-mean-square value was calculated for different values of the number of vortex structures, N v , in the computational domain, with fixed value of N V Γ 2 = 2000. The standard deviation and mean values of these results were obtained, the ratio of which yields the relative standard deviation σ E /E and σ Ω /Ω. A plot of the relative standard deviations is shown in Figure 3(b) as functions of N V . The standard deviations exhibit some variation with the number of vortices for small values of N V , but for N V ≥ 64 they are nearly independent of the number of vortices. The standard deviation for turbulent kinetic energy is about 4%-5% of the mean value, whereas that for enstrophy is only about 1% of the mean value. We note that this deviation is not a resolution error; since enstrophy is computed from the velocity gradients it is significantly more sensitive to resolution errors than is the kinetic energy field. Rather, the observed fluctuations arise from the variation in position and orientation of the vortices between the different configurations examined. Since the vorticity field is largely confined to the region within and immediately surrounding the vortex structures, it is reasonable that the relative standard deviation for enstrophy should be small, provided that the vortex structures do not overlap. The higher value of the relative standard deviation for turbulent kinetic energy arises from the fact that the velocity field at any point in the flow is dependent not only on its position relative to the nearest vortex structure but rather on all vortex structures in the flow field.
The power spectrum e(k) was examined for a series of computations with N V Γ 2 = 200 and numbers of vortices of N V = 512, 256, 128, 64, and 32, with values of Γ adjusted to give the specified product value. The spectrum lines fall on top of each other and cannot be distinguished, which confirms the prediction from (25) that the spectrum depends on N V and Γ through the combination N V Γ 2 ℓ/V . In Figure 4, we compare the SVS computational spectrum for the case with N V = 512 to Saffman's approximate prediction (25). The theoretical expression is found to be significantly higher than the SVS computational values, particularly for higher values of k. This result is likely due to the fact that Saffman assumed all vortices to be non-interacting, and so he simply added the kinetic energy of each vortex (associated with its own self-induced velocity) to obtain the total kinetic energy. In the computations, the vortex structure orientation is random, so the induced velocity from one structure will counter that from other structures at sufficiently large distances, thus reducing the total kinetic energy. Also shown in Figure 4 is a line indicating k −5/3 dependence, which fits the computational plot reasonably well within the low-wavenumber inertial range, similar to the observations of Kivotides and Leonard (2003).
The SVS predictions are compared in Figure 4 to the results of a pseudo-spectral direct numerical simulation (DNS) computed on a 128 3 grid, similar to that presented by Vincent and Meneguzzi (1991). The flow is initiated by a randomly perturbed velocity field with uniform probability distribution for wavenumbers spanning the interval 1 ≤ k ≤ 64. Dealiasing is performed by setting the coefficients of the highest 1/3 wavenumber coefficients to zero using a spherical filter. A preliminary computation is run without forcing until time t = 10 in order to allow the turbulence to develop a range of length scales characteristic of statistically stationary homogeneous isotropic turbulence. The computation is then restarted with non-zero forcing, where the transform of the forcing vector is assumed to be proportional to the fluid velocity transform, such that (Lundgren, 2003 andRosales andMeneveau, 2005)  where the coefficient C is adjusted at each time step so as to maintain approximately constant turbulent kinetic energy. The current computations are performed with k crit = 5, so that the forcing acts only on the large-scale eddies. Various parameter values characterizing the DNS computations are given in Table II. The spectrum predicted by the DNS computations compares well with the SVS predictions for low values of wavenumber (k < 20), but for high wavenumber the SVS spectrum decays much more quickly than do the DNS results. This rapid decay at high wavenumber is consistent with the fact that the vortex radius for these computations was specified to be eight times larger than the Kolmogorov length scale, so the SVS flow field has little energy at very small length scales. The velocity probability density function (PDF) in one coordinate direction (x-direction), normalized by the root-mean-square value, was computed for a series of SVS computations with N V Γ 2 = 3975 and different number of vortex structures. Unlike the power spectrum, the velocity PDF exhibits significant variation with a value of N V . This observation indicates that the velocity PDF varies with N V and Γ independently, and not only through the product N V Γ 2 . The PDF has a fat tail for low values of N V , typical of a superstatistical system (Beck, 2008), but the PDF functions for large values of N V (greater than about 500) approach an asymptotic curve that is nearly Gaussian. In Figure 5(a), a comparison is shown of the velocity PDF for the SVS computation with N V = 512, the DNS simulation (symbols), and a best-fit Gaussian curve p(v) = 0.8 exp(−0.5v 2 ), where v ≡ v x /v x,r ms . The DNS results are in close agreement with the Gaussian function, as expected (Voth et al., 1998). The SVS predictions fit well to the Gaussian function for v x /v x,rms < 3, but for higher values of v x they exhibit higher values. This difference indicates that while still very rare, high velocity occurrences are more common for the SVS computations than for the DNS.
The PDF of the x-component of the fluctuating fluid acceleration field is plotted in Figure 5(b). Fluid acceleration is computed from the SVS or DNS velocity field for post-processing purposes using a centered difference approximation in space and a forward difference in time. We again find that the PDF plot is sensitive to the value of N V , and that it approaches an asymptotic curve for values of N V greater than about 500. The SVS prediction for the case with N V = 512 is compared to the DNS results in Figure 5(b). Also shown in this figure is the empirical expression for the PDF, obtained experimentally by La Porta et al. (2001). In this expression, a ≡ a x /a x,r ms , and the coefficients are given by a best fit to the experimental data of La Porta et al. as c 1 = 0.539, c 2 = 1.588, and c 3 = 0.508. The SVS prediction for acceleration PDF with N V = 512 is found to agree closely with both the DNS prediction and with the experimental correlation (28), and in all cases the acceleration PDF exhibits non-Gaussian statistics characterized by fat tails, typical of a highly intermittent signal. Mordant et al. (2004) associates the acceleration intermittency in turbulent flows with the presence of coherent vortex structures, so agreement between the SVS and DNS for the acceleration PDF is another indication that the coherent vortices are correctly modeled in the SVS representation. The PDF of the vorticity component ω x is plotted in Figure 6 from SVS results with N v = 2048 vortex structures in the computational domain. The vorticity is determined by first computing the synthetic turbulence velocity field, as discussed in Section III, and then numerically differentiating using a centered finite-difference method to obtain vorticity from ω = ∇ × u. The PDF for vorticity is sensitive to the number of vortex structures used for the SVS computations, and because the vorticity is evaluated using a velocity gradient it required a somewhat larger number of vortices to reach the asymptotic state for large vortex numbers than did the velocity or acceleration PDFs. The SVS vorticity PDF is shown in Figure 6 to be in excellent agreement with the vorticity PDF obtained from the DNS predictions.

V. VALIDATION OF SVS PREDICTIONS FOR PARTICLE COLLISION RATE
The SVS predictions for particle interactions were validated by comparison to DNS results with use of a soft-sphere discrete-element method (DEM) simulation for a set of N p colliding non-adhesive particles of diameter d and mass m. The computations solve the momentum and angular momentum equations for the particle velocity and rotation rate, given by subject to forces and torques induced by the fluid flow (F F and M F ) and by the particle collision and adhesion (F A and M A ). Here, I is the moment of inertia and v and Ω are the particle velocity and rotation rate, respectively. The dominant fluid force is the particle drag force, but we also accounted for secondary forces including the Saffman and Magnus lift forces and the added mass and pressure gradient force on the particles. Particle Reynolds numbers were small, allowing use of the Stokes drag law and low Reynolds number lift laws (Saffman, 1965(Saffman, , 1968and Rubinow and Keller, 1961). Collisions were detected when the distance between two particles is less than the particle diameter. Collision forces between the particles include the normal elastic and dissipative forces, sliding resistance, and twisting resistance. Particle normal collision was computed for non-adhesive particles using the nonlinear Hertz (1882) theory for normal elastic force, the Tsuji et al. (1992) model for normal dissipative force, and the Cundall and Strack (1979) model for sliding resistance. The fluid velocity was interpolated from a 128 3 fluid grid onto the Lagrangian particle locations with cubic accuracy using the M4 ′ variation of the B-spline interpolation method, which was originally developed by Monaghan (1985a) and is commonly used in spherical particle hydrodynamics (Monaghan, 1985b) and for regridding in vortex methods (Cottet and Koumoutsakos, 2000). The multiple time step algorithm of Marshall (2009) was used with three different time step levels, corresponding to the fluid, particle, and collision time scales, arranged from largest to smallest. The reported computations used a fluid time step of ∆t = 0.01 for a duration of 10 000 time steps. The DNS runs were initiated using a preliminary computation without particles with 5000 time steps to establish a statistically-steady turbulent flow. Simulations were performed on a cubic grid with 2π side length and 46 656 particles. A listing of integral flow measures for the different cases examined in this comparison is given in Table III. The number of vortices was varied from N V = 32 to 2048 in the SVS runs SVS-1a through SVS-1g in order to examine the effect of the number of vortices on the collision results, and in each case the value of vortex circulation was adjusted to maintain nearly constant turbulent kinetic energy. The computations were performed for values of the integral-scale Stokes number St 0 of 0.07, 0.34, and 1.7, where St 0 is defined by (1) with u ℓ = u 0 and ℓ = ℓ 0 . The corresponding values of the Kolmogorov-scale Stokes number St K for these three cases are 0.81, 3.94, and 19.9, respectively. A filtered DNS computation (DNS-F) was also performed in which the coefficients of the highest 67% of the wavenumbers (k > 21.3) was set to zero, which yields an energy spectrum very close to the SVS spectrum. The filtered DNS run is used as a method to determine the influence of small-scale fluctuations on the particle collisions. Besides kinetic energy, integral measures listed in Table III include enstrophy Ω, a vorticity magnitude measure ω 95 , and a stretching rate measure S. The vorticity magnitude measure ω 95 is defined as the value of vorticity magnitude for which 95% of the grid points have a lower vorticity magnitude. The stretching measure S is defined as the average over the flow field of the maximum value of the logarithmic stretching rate λ 1 =Λ/Λ. Here, Λ is the stretch of a material line segment along the principal direction of the rate of deformation tensor D associated with the largest eigenvalue λ 1 of D. Since D is symmetric, the eigenvalues of D can be efficiently computed using the Smith algorithm (Smith, 1961). The enstrophy for the filtered DNS run (DNS-F1) is about twice the value for the associated SVS run (SVS-1), and the enstrophy for the unfiltered DNS run (DNS-1) is about 20% higher than that for the filtered DNS run due to the contribution of the small vortices filtered out in the DNS-F1 run. In accordance with the result (24), the enstrophy remains nearly constant in the SVS runs (SVS-1a through SVS-1g) as the number of vortices is changed with N V Γ 2 held constant. The vorticity magnitude parameter ω 95 is about 40% larger and the stretching measure S is about 15% larger for the DNS run compared to the SVS-1a run.
The total number of collisions was found to increase almost linearly with time, and the slope of this line was used to compute the collision rate per unit volumeṅ C . From this value, the collision kernel α 11 was computed using the definitioṅ where n = N p /V is the number of particles per unit volume. The predicted value of α 11 for each case was computed from (30) using the specified value of n and the computed value ofṅ C based on a linear fit to the total number of collisions, and the resulting values of collision kernel are listed in Table III. A comparison of the collision kernels between the full DNS, the filtered DNS, and the SVS method was conducted for integral-scale Stokes numbers of St 0 = 0.07, 0.34, and 1.7, where the Stokes number is changed by modification of the particle diameter. As predicted by collision theory Turner, 1956 andAbrahamson, 1975), the collision kernel increases with particle diameter (indicated by increasing Stokes number), with DNS predictions of α 11 = 5.8 × 10 −5 , 3.26 × 10 −4 and 6.45 × 10 −3 for St 0 = 0.07, 0.34, and 1.7, respectively. The filtered DNS predictions for the collision kernel are within about 4% of the full DNS predictions for each case, indicating that the small scales of the turbulent motion have little effect on the collision coefficient. The collision kernel for the SVS model with 2048 vortices was about 16% lower than the full DNS prediction for the St 0 = 0.34 case, and the SVS model predictions for α 11 were within 0.8% and 5.7% of the full DNS predictions for the St 0 = 0.07 and 1.7 cases, respectively. The effect of the number of vortex structures on the SVS predictions was examined by repeating the run for St 0 = 0.34 with 32, 64, 128, 256, 512, 1024, and 2048 vortices, while at the same time adjusting the vortex strength to keep the kinetic energy approximately constant. The tendency of particles to cluster can be characterized by the radial distribution function (RDF), g(r), defined by where the average number of particles per unit volume ρ 0 is related to the particle volume fraction C p by ρ 0 = 6C p /π, and N(r) is the average number of neighboring particles whose centroids are located within a radial distance r from a given particle centroid. The value of g(r) is estimated by counting for each particle the number of neighboring particles that fall into a set 400 spherical bins, each of width ∆r = 0.000 15, surrounding the given particle. The number of particles in each bin is averaged over all particles in the computational domain and over 1000 time steps near the end of the computations in order to smooth the distribution. Figure 7 shows a comparison of the RDF for both As noted by Zaichik et al. (2006), the collision kernel is proportional to the product of the radial distribution function g(r) (RDF) and the relative radial velocity ⟨w r ⟩ (RRV) evaluated at collision (r = 2r p ). Each of these quantities was separately computed for cases with different Stokes number to examine the individual quantities that make up the collision kernel. A set of plots is given in Figure 8 showing RDF and RRV at collision as a function of the Kolmogorov-scale Stokes number St K for both DNS results and SVS results with N V = 2048 vortices. Our predictions are compared to the DNS results for RDF and RRV of Fayed and Ragab (2013) for Re λ = 77 and of Wang et al. (2000) for Re λ = 75, and to the DNS results for RDF of Sundaram and Collins (1997) for Re λ = 54. The RDF value for St K = 19.9 is nearly the same in the SVS and DNS predictions, and so the two symbols for DNS and SVS results are almost coincident in Figure 8(a). The RDF exhibits a very thin peak near the collision point for small Stokes number, which contributes to the high data variability in Figure 8 Wang et al. (2000) at Re λ = 75 (circles, red line), DNS data of Sundaram and Collins (1997) for Re λ = 54 (squares, blue line), DNS data of Fayed and Ragab (2013) for Re λ = 77 (deltas, green line), and our DNS predictions (filled diamonds) and SVS predictions (open diamonds) for Re λ = 81.

VI. VALIDATION OF SVS PREDICTIONS FOR TURBULENT PARTICLE AGGLOMERATION
Computations to examine turbulent agglomeration were conducted with a similar DEM algorithm as described in Sec. V, but with a modification of the collision force and torque models to account for adhesion effects. In particular, the normal elastic and adhesive van der Waals force was computed using the model of Johnson et al. (1971) (i.e., the JKR model). Adhesion introduces a strong rolling resistance torque, for which we used the model of Dominik and Tielens (1995), along with experimental results of Ding et al. (2008) to set the critical angle for onset of particle rolling. The effect of adhesion on the sliding resistance was modeled using an expression derived by Thornton (1991). We also included a crowding correction term for the particle drag force developed by Di Felice (1994). A comprehensive summary of the computational method for both adhesive and non-adhesive particles is given by Marshall (2009). The reported computations used a fluid time step of ∆t = 0.005 for a duration of 20 000 time steps, with a total of 46 656 particles. As discussed in Sec. V, the DNS runs were initiated using a preliminary computation without particles with 5000 time steps to establish a statistically-steady turbulent flow.
A particle agglomerate constitutes a set of particles which are bonded to each other, either directly or via other intermediate particles of the agglomerate, via soft (e.g., van der Waals) bonds. A set of particles bonded via hard bonds (e.g., sintered particles) is referred to as an aggregate, and is outside the scope of this paper. Agglomerate development in the turbulent flow field is characterized in the current paper using two dimensionless parameters-the Stokes number St and the adhesion parameter Ad. The adhesion parameter Ad is a ratio of adhesive force to particle inertia, defined by the ratio where the adhesive surface energy density γ is equal to half the work required to separate two surfaces that are adhesively bound per unit surface area. Both the Stokes number in (1) and adhesion parameter were defined using the characteristic length scale ℓ 0 and velocity scale u 0 of the turbulence integral scale for the fluid length and velocity scales ℓ and u ℓ , which is indicated by a subscript "0." Plots showing SVS predictions for the total number of particles contained in an agglomerate, N tot , and the average number of particles per agglomerate, N pagg , as functions of time are given for a case with St 0 = 0.34 and Ad 0 = 11 in Figure 9 for different values of the number of vortex structures, N V , ranging from 128 to 2048. The vortex strength is adjusted to maintain a constant turbulent kinetic energy in each case. While the collision kernel listed in Table III  The run shown in Figure 9 was extended to a time of t = 250 to examine the continued agreement between SVS and DNS as the equilibrium condition is reached. The average number of particles per agglomerate is plotted versus time for this extended run in Figure 10, showing that the SVS run (with N V = 2048) and DNS continue to exhibit reasonable agreement at long time. The value of the average number of particles per agglomerate fluctuates in time when this statistical equilibrium state is reached due to breakup and recollision of large agglomerates.
The effect of Stokes number is shown in Figure 11, which compares SVS predictions with N V = 2048 vortex structures and DNS predictions for values of the Stokes number of St 0 = 0.1, 0.2, and 0.34, in all cases with Ad 0 = 11. The different Stokes numbers are produced by changing the particle diameter, with all other parameters held constant. Plots are given both for the average number of particles per agglomerate, N pagg , and for the total number of agglomerates, N agg , as functions of time. The value of N pagg decreases rapidly with decrease in St 0 , going from N pagg = 130 at t = 100 for St 0 = 0.34 to N pagg = 10 for St 0 = 0.1. The peak value of the number of agglomerates is shown in Figure 11(b) to be nearly the same for the three cases, but the peak occurs at a later time as the Stokes number decreases. The observed differences in agglomeration measures with change in St 0 are primarily due to decreasing collision rate as the Stokes number decreases, which is consistent with theoretical predictions for collision rate at both small and large Stokes numbers Turner, 1956 andAbrahamson, 1975). Good agreement is observed between the SVS and DNS predictions.
The effect of adhesion parameter is examined in Figure 12, which compares SVS predictions with N V = 2048 vortex structures and DNS predictions for values of the adhesion parameter of Ad 0 = 5.5, 11, 28, and 110, in all cases with St 0 = 0.34. The different adhesion parameter values are produced by changing the adhesion surface energy density γ, with all other parameters held constant. As expected, the average number of particles per agglomerate decreases in Figure 12(a) with decrease in Ad 0 . The total number of agglomerates in Figure 12(b) is found to peak at nearly the same time for the different values of Ad 0 , and to then decrease rapidly after the peak value for high values of Ad 0 , indicating that agglomerates are colliding to form larger agglomerates. For Ad 0 = 5.5, the number of agglomerates decreases slowly after the peak since colliding agglomerates might not adhere to each other or might breakup again into smaller agglomerates. Again, good agreement is observed between the SVS and DNS predictions.
The agglomerate number distribution indicates the percentage of agglomerating particles contained in agglomerates consisting of n particles. The agglomerate number distribution is sorted into logarithmic bins of base 2, where the value of bin size indicates the nominal number of particles in  Figure 13. SVS predictions with N V = 2048 vortex structures are observed to yield a number distribution that is reasonably close to that obtained using DNS.
Each agglomerate is characterized by the number of particles N contained in the agglomerate and the radius of gyration R g , which is defined by In this equation,x denotes the position vector of the agglomerate centroid and x i is the centroid of the ith particle within the agglomerate. It is well known that particle agglomerates admit a power law relating N and R g given by (Adachi and Ooi, 1990;Liu et al., 1990;and Jiang and Logan, 1991) where K is a coefficient (called the fractal pre-factor) and the exponent D is called the fractal dimension of the agglomerate. The value of D varies over the interval 1 ≤ D ≤ 3 depending on the agglomeration formation mechanism (Brasil et al., 2001). For instance, Eggersdorfer et al. (2011) cite typical values of D = 2.5 for diffusion-limited agglomeration, D = 3.0 for ballistic particle-cluster agglomeration, and D = 1.8 for diffusion-limited cluster-cluster agglomeration. For turbulent agglomeration of latex particles in stirred tanks, Selomulya et al. (2001) report values of D between 1.7 and 2.1 and Waldner et al. (2005) report values of D between 1.8 and 2.6. A log-log plot of N versus R g /d for both DNS results and SVS predictions with N V = 2048 is shown in Figure 14. The DNS and SVS predictions are in excellent agreement, and both are found to exhibit a best-fit line with slope D = 2.3. As discussed above, this value of fractal dimension of the FIG. 14. Plot showing the number of particles in an agglomerate N versus the ratio of the gyration radius to the primitive particle diameter, R g /d, with both DNS data (triangles, blue) and SVS data with N V = 2048 (crosses, red). The solid line is a best-fit to the data with a slope of D = 2.3. particle agglomerates is in good agreement with values noted in previous experimental literature for turbulent agglomeration.

VII. CONCLUSIONS
An accelerated form of the stochastic vortex structure (SVS) method for subgrid-scale turbulence modeling for interacting particles was developed using the method of multipole expansions. It was shown that with only five box levels, the accelerated method can reduce the velocity computation time by two orders of magnitude, with error in the total turbulent kinetic energy (TKE) prediction of less than 2%. The effect of the stochastic nature of the SVS algorithm on prediction of mean quantities was examined, and it was found that the ratio of the standard deviation to the mean value obtained from repeated runs with different vortex positions and orientations was about 5% for TKE and 1% for enstrophy. Characteristics of the SVS synthetic turbulence predictions were examined against results of direct numerical simulation (DNS) and various theoretical and experimental results from the literature. The predicted energy spectrum was compared against both DNS results and approximate theoretical results from Saffman (1997), and shown to be in reasonable agreement with both for moderate and small values of wavenumber (less than about 20), but (as expected) to give too low values for higher wavenumbers. The predicted velocity, acceleration, and vorticity probability density functions (PDFs) were found to be sensitive to the number of vortex structures used, but to approach the DNS predictions for large number of vortex structures. SVS predictions for other integral measures, such as the ω 95 measure of the maximum vorticity magnitude and the average stretching rate measure, also exhibit good agreement with DNS.
Though the validation of the statistical properties of SVS-generated fields is encouraging, the ultimate arbiter of the robustness of this model is whether or not it achieves the ultimate modeling objectives. To this end, simulations with colliding, non-adhesive particles were performed comparing the SVS predictions for radial distribution function, relative radial velocity, and collision kernel to DNS results. Computations were performed for an integral scale Stokes number range of 0.07-1.7, yielding good agreement between SVS and DNS predictions. The simulations indicate that the SVS results for collision rate are not very sensitive to the number of vortex structures as long as this number is sufficiently large. DNS and SVS simulations were also performed for collision and agglomeration of adhesive particles over a range of Stokes number and adhesion parameter values. Agglomeration measures examined include total number of particles captured in agglomerates, number of agglomerates, average number of particles per agglomerate, number distribution of agglomerates, and agglomerate fractal dimension. Values of these agglomeration measures were found to approach values close to those of the DNS predictions for sufficiently high numbers of vortex structures.
The paper suggests that the stochastic vortex structure method provides a rapid, reliable approach for modeling subgrid-scale turbulence fluctuations for flows with interacting particles. The SVS method is consistent with the large-scale energy spectrum and the various probability density function curves that describe homogeneous turbulence, as well as with a wide range of integral measures of the turbulent flow. The speed-up in velocity field computation introduced in the current version of the SVS method makes this approach highly efficient compared to other synthetic turbulence approaches. Because the SVS method deals directly with the vortical structures that dominate the large-scale motion of the turbulence, it allows accurate prediction of phenomena, such as particle clustering, that are dependent on the structural form of the turbulent eddies.
We note that the current validation study was conducted for a relatively low Reynolds number flow for which the integral-scale Stokes number was close to unity. For high Reynolds number turbulence, there exists a large range of scales between the integral scale and Kolmogorov scale. A study using the wavelet-based coherent vortex simulation approach by Nejadmalayeri et al. (2013) found that the number of energy-containing structures at a fixed kinetic energy level increases linearly with Reynolds number in homogeneous turbulence. While the SVS method has not yet been tested for high Reynolds numbers for purposes such as prediction of particle collision rate, we speculate that it may not be necessary to cover the entire range of these length scales with the synthetic turbulent flow. Rather, it might be sufficient to introduce SVS structures only for a length scale ℓ for which the eddy Stokes number St ℓ is closest to unity. Eddy structures much larger than this scale ℓ will simply advect the particles with minimal relative motion between the particles, and the fluctuations associated with eddies much smaller than ℓ will be filtered out by the particle inertia. However, we also recall that several experimental and computational studies have observed that intense vortex structures are less prominent for high Reynolds number turbulence (with Re λ ≥ O(1000)) than is the case at low Reynolds numbers (Berlin et al., 1996 andIshihara et al., 2009). The potential effectiveness of vortex-based methods such as SVS at high turbulent Reynolds numbers will therefore need to be carefully assessed in future work.

ACKNOWLEDGMENTS
This research was supported by the U.S. National Science Foundation under Grant No. CBET-1332472.

APPENDIX: INDUCED VELOCITY BY A SCALAR GRADIENT VORTICITY FIELD
Substituting (4) into the Biot-Savart Equation (6) gives the induced velocity at a point x as Making use of the identity s/s 3 = −∇(1/s) = ∇ ′ (1/s) and the vector identity ∇ × ∇ζ = 0, Green's theorem can be used to write the integral in the second term on the right-hand side of (A1) as 1 4π where S is the bounding surface of V . At large distances, |x| >> L, the gradient field ∇ζ has the form of a dipole that decays with distance r as O(1/r 3 ). Consequently, the surface integral in (A2) approaches zero as S → ∞, leading to the conclusion that the velocity field is entirely induced by the non-gradient part ω * of the vorticity field.