A stochastic vortex structure method for interacting particles in turbulent shear flows

In a recent study, we have proposed a new synthetic turbulence method based on stochastic vortex structures (SVSs), and we have demonstrated that this method can accurately predict particle transport, collision, and agglomeration in homogeneous, isotropic turbulence in comparison to direct numerical simulation results. The current paper extends the SVS method to non-homogeneous, anisotropic turbulence. The key element of this extension is a new inversion procedure, by which the vortex initial orientation can be set so as to generate a prescribed Reynolds stress field. After validating this inversion procedure for simple problems, we apply the SVS method to the problem of interacting particle transport by a turbulent planar jet. Measures of the turbulent flow and of particle dispersion, clustering, and collision obtained by the new SVS simulations are shown to compare well with direct numerical simulation results. The influence of different numerical parameters, such as number of vortices and vortex lifeti...


I. INTRODUCTION
Computational modeling of the motion of interacting particles, droplets, or bubbles subject to subgrid-scale fluctuations in turbulent flows is a long-standing challenge in multiphase flow simulations. The Reynolds-averaged Navier-Stokes (RANS) approach remains the most common method for engineering solution of practical turbulent flows, providing both manageable computation times and reasonably accurate prediction of key flow features, such as boundary layer separation. However, when used in conjunction with the Lagrangian simulation of particulate fluids, it is necessary to augment the RANS equations with some model to account for the effect of the turbulent fluctuations when computing the particle trajectories. This problem also arises when using the large eddy simulation (LES) approach with sufficiently small values of the Kolmogorov-scale particle Stokes number. 1 A number of effective methods exist for dealing with this problem for non-interacting particles (see Ref. 2 for a review), but subgrid-scale modeling for transport of interacting particles in turbulent flows remains an unresolved challenge. Particle interaction is essential in a wide range of turbulent flow problems occurring in nature, such as turbulence effects on collision of rain droplets or snow flakes, contact electrification of dust particles in sand storms, and agglomeration of particles in volcanic plumes or of pollution particulates in the atmosphere. Particle interaction also plays an important role in many industrial particulate flow problems, such as pharmaceutical manufacturing, paint production, wastewater treatment, additive manufacturing processes, 3D printing, flame synthesis of nanoparticles, and fly ash capture from combustion furnaces. a) Author to whom correspondence should be addressed: jmarsha1@uvm.edu.
The challenges associated with subgrid-scale modeling for turbulent transport of interacting particles arise from three considerations. First, it is critical for simulation of particle interaction to accurately model small values of the particle separation distance. However, small separation distances imply that the nearby particles are experiencing forcing from the same set of nearby turbulent eddies so that the fluid fluctuation velocity at the particle positions is highly correlated. As a consequence, any model in which each particle experiences uncorrelated forcing will not be appropriate for simulation of interacting particles. Second, particle collision and adhesion processes occur over time scales that are very small, typically much smaller than those associated with the fluid flow. The numerical calculation consequently becomes numerically stiff when particle interactions are included, particularly when using methods such as the soft-sphere discrete-element method (DEM). Synthetic turbulence models commonly used to approximate the subgrid-scale turbulent fluctuations must therefore be highly efficient in order to be manageable with small time steps. Third, turbulent eddy structures are known to expel particles with higher density than the surrounding fluid, leading to formation of particle clusters in the region in between the eddies that can have local particle concentrations an order of magnitude or more above the average concentration. [3][4][5][6][7] This phenomenon leads to the so-called preferential concentration effect, which can dramatically increase the particle collision rate, agglomeration, and other interactions in these high-concentration regions. 8,9 Since particle clustering in turbulent flows occurs due to interaction of particles with coherent eddies, it is natural to utilize a vortex structural approach in modeling the effect of turbulent fluctuations on interacting particles. Vortex structural models have long been used in turbulence flow modeling, dating back to Townsend's 10 model of homogeneous turbulence as a collection of Burger's vortices and Lundgren's 11 spiral vortex model, as well as extensions of these models by Pullin and Saffman 12 and Saffman. 13 The scaling and structure of coherent vortices were studied numerically for homogeneous turbulence by Jiménez et al. 14 and experimentally for turbulent shear flows by Berlin et al., 15 among others. Both studies found that the turbulent vorticity field 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) vorticity. The length and core radius of the coherent vortices were found to scale with the Lagrangian 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. Theoretical proof of these scaling observations was provided by Kambe and Hatakeyama. 16 Kivotides and Leonard 17 reported a computational study in which homogeneous turbulence is represented by a set of finite-length vortex structures and showed that this system generates an energy spectrum that satisfies the Kolmogorov k −5/3 scaling in the turbulence inertial range. The effectiveness of vortex structural models for prediction of turbulence structure functions and various velocity and vorticity probability density functions was discussed by Refs. 18-22. Extensions of the vortex filament method were successfully utilized for simulation of a number of turbulent shear flows, including mixing layers, 23 co-flowing jets, 24 and boundary layers. 25 Ayyalasomayajula et al. 26 proposed a vortex structural model for transport of particles in homogeneous isotropic turbulence using a two-dimensional array of uniformly spaced vortices, where a stochastic algorithm is used to determine the vortex strength. Somewhat surprisingly, given the highly simplified nature of this model, the predicted particle acceleration statistics and clustering were similar to direct numerical simulation (DNS) results. Sala and Marshall 27 proposed a three-dimensional stochastic vortex structure (SVS) model, again for homogeneous isotropic turbulence, where the turbulent eddies are represented by a set of finite-length vortex structures which are randomly positioned and oriented in the flow field. The vortex length and core radius were assumed to be proportional to the turbulence integral and Kolmogorov length scales, respectively. Unlike the vortex filament method, the SVS method does not use the vortex structures to evolve the turbulent flow field; instead, the vortex structures are used only to approximate a subgrid-scale synthetic turbulence to use for particle evolution in a flow with a given Reynolds stress distribution. An accelerated version of the SVS method was developed by Dizaji and Marshall 28 using both the fast multipole method and a local Taylor series expansion which speeds up the computations by up to two orders of magnitude with negligible difference in flow field or particle interaction statistics. The SVS model was shown to yield predictions for turbulence energy spectrum, velocity and acceleration probability density function (PDF), and particle collision rate that are in close agreement with DNS predictions. Dizaji and Marshall 28 also verified that the SVS model is highly effective at accurately predicting various measures characterizing agglomerate formation for adhesive particles in turbulent flows.
One criticism of the SVS model is that, to date, all applications of this model have been for isotropic, homogeneous turbulence. The objective of the current paper is to extend the SVS model to non-homogeneous, anisotropic turbulent flows and to validate this extended model in comparison to direct numerical simulation (DNS) results. The extension of the SVS model for anisotropic turbulence is described in Sec. II A, with particular focus on a proposed inversion algorithm by which the orientation of the SVS vortex structures can be adjusted to yield a prescribed Reynolds stress field. This vortex structure initialization method is examined and validated in Sec. II B for both homogeneous and inhomogeneous anisotropic flow fields. Computational methods used for particle transport and for direct numerical simulation (DNS) for validation of the SVS model are described in Sec. III. Validation of the SVS model with comparison to DNS results for particulate turbulent planar jet flow is reported in Sec. IV. Conclusions are given in Sec. V.

A. Anisotropic SVS method
The stochastic vortex structure (SVS) model approximates the turbulent vorticity field by a collection of vortex structures placed in the flow field. In its simplest version, the vortex structures in the SVS model all have the same finite length L, core radius δ, and strength Γ. The vortex length L is of the order of magnitude of the turbulence Lagrangian integral length scale 0 = 0.5 u 0 3 /ε, where u 0 is the turbulence rootmean-square velocity and ε is the turbulence dissipation rate per unit mass. The core radius δ of the coherent vortices was estimated numerically by Jiménez et al., 14 experimentally by Berlin et al., 15 and theoretically by Kambe and Hatakeyama 16 to be 3-4 times the Kolmogorov length scale, η = (ν 3 /ε) 1/4 , where ν is the fluid kinematic viscosity. In the current work, we use somewhat larger vortex structures with a core radius of 8η so as to ensure sufficient number of grid points to resolve the velocity variation across the vortex cores; however, SVS computations were repeated with core radius δ = 4η, and the results were found to be almost identical to those with larger core radius. Each vortex structure has a lifetime T V which is proportional to the integral time scale, T = q/3ε, where q = 1.5 u 2 0 is the turbulent kinetic energy per unit mass, although we note that the model results are not sensitive to the choice of vortex lifetime.

Vortex structure initialization
The Reynolds stress tensor R has components in the global Cartesian coordinate system given by R ij = u i u j , where a prime denotes the fluctuating velocity component and an overbar denotes a time average. In the SVS simulation, the anisotropy of the turbulent fluctuations is produced via preferential orientation of the vortex structures. It is necessary to develop a method for specifying the probability distribution of the vortex structure orientation so as to be consistent with the given Reynolds stress tensor, which is a type of inverse problem. Turbulence anisotropy is related both to differences in value of the three normal components of the Reynolds stress and to the offdiagonal Reynolds stress components. We employ a four-step approach for setting the vortex orientation in accordance with a given Reynolds stress tensor, as described below. Prior to implementing this procedure, we compute a set of M = 642 evenly spaced test points on the surface of a unit sphere by dividing the faces of an icosahedron a prescribed number of times and projecting the vertices to the unit sphere.
In the first step, the Reynolds stress tensor is interpolated from the grid covering the flow field onto the centroid position of a vortex structure. In the second step, we rotate the coordinate system to a frame in which the Reynolds stress tensor at the vortex centroid is diagonal. This is achieved by computing the set of three eigenvalues λ (k) and associated normalized eigenvectors x (k) of the Reynolds stress tensor. We define a principal direction coordinate system as a coordinate frame whose base vectors are the three eigenvectors of R. The components of the Reynolds tensor in the principal direction coordinate system, denoted by R * ij , are given by . (1) In the third step, the vortex structure orientation is set in this principal direction coordinate system. The inverse procedure by which this is achieved is based on the observation that a vortex structure oriented in the x-direction, say, would induce a velocity field in which R 11 = 0 and R 22 = R 33 . We define vortex orientation weighting coefficients c 1 , c 2 , and c 3 , normalized by c 1 + c 2 + c 3 = 1, such that Solving system (2) for the three orientation weighting coefficients gives The orientation of a vortex structure is specified at the initial time step by randomly selecting one of the M test points on the unit sphere, obtained using the procedure described at the beginning of this section. The coordinates of the selected test point in principal direction coordinates are denoted by (ξ * 1 , ξ * 2 , ξ * 3 ). Using the weighting coefficients obtained in (3), the vortex structure orientation ζ is set in principal direction coordinates as where ∆ ≡ [(c 1 ξ * 1 ) 2 + (c 2 ξ * 2 ) 2 + (c 3 ξ * 3 ) 2 ] 1/2 . In isotropic turbulence, the three orientation weighting coefficients are equal so that (4) results in random vortex orientation with uniform distribution.
The fourth step of the vortex structure initialization process is to rotate the structure back into the global coordinate system used for the computation. We recall that the components of the rotation tensor A from the global Cartesian coordinates to a principal direction coordinate system form an orthonormal 3 × 3 matrix whose three columns are the components of the three eigenvectors x (k) . The components of the vortex structure orientation vector ζ in the global coordinate frame can therefore be written in terms of the components in (4) as

Velocity calculation
The vortex structures induce a velocity field u, which is computed on the flow grid using the fast multipole acceleration method described by Dizaji and Marshall. 28 The accelerated velocity computation method first partitions the computational domain into a tree-structure composed of uniform-size boxes, where at every level of the tree structure each box from the previous level is divided into eight "offspring" boxes by dividing the side lengths in half in each direction. For each of the smallest "target" boxes in the tree structure, a list of other "source" boxes with which it interacts "directly" and "indirectly" is developed based on the analytical error estimate for the multipole expansion by Salmon and Warren. 29 For source boxes on the direct interaction list, we compute the induced velocity from each vortex structure in the source box on each grid cell node in the target box by interpolation from a planar section, where the induced velocity from a unit strength vortex structure on the plane is pre-computed at the start of the simulation. For source boxes on the indirect interaction list, the induced velocity from all sources in the box is computed at the centroid of the target box using multipole expansion. 30,31 The contribution of this induced velocity at the individual grid cell nodes within the target box is then determined using a local Taylor series expansion. Induced velocity from vortex structures from one period of the computational domain in each direction is also included in the computation. This accelerated method was shown by Dizaji and Marshall 28 to produce very accurate results with computation times that are nearly two orders of magnitude less than the direct computation method using only four levels of the box structure.

Vortex time evolution
Each of the N V vortex structures are 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. The initial age of the n th 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 x n and orientation given by the same four-part procedure as used to initialize the vortex structure orientation.

Limitations of inversion method
The inverse method for initialization of the SVS vortex structures described in Sec. II A is validated in this section for different test computations in which the vortex structures are initialized using a prescribed Reynolds stress field, and then the Reynolds stress is evaluated from the computed SVS velocity field and compared to the prescribed field. In conducting this validation, it is important to bear in mind that the inversion procedure described in Sec. II A is subject to limitations, which can be ultimately associated with the fact that we are attempting to generate a turbulence field using only the induced velocity from tubular vortex structures. Mathematically, these restrictions require that the three coefficients c 1 , c 2 , and c 3 defined in (3a)-(3c) must all be positive. This in turn introduces the following three restrictions on the values of the eigenvalues λ (k) : If we now consider the special case of a two-dimensional turbulent mean flow, such as a plane jet or channel flow, the Reynolds stress tensor R ij has the form Solving for the eigenvalues of the Reynolds stress tensor gives Using R 33 for normalization, Reynolds stress ratios can be defined as which are all positive by definition. From the solutions (10), we find that the limitation (8b) is always satisfied and the limitations (8a) and (8c) become, respectively,

Special case Prescribed values Limitation
Specific limitations for several special cases, as computed from (12), are listed in Table I.

Validation for homogeneous turbulence
The inversion method described in Sec. II A was validated first for the case of homogeneous turbulence, in which the Reynolds stress is uniform in space. The Reynolds stress tensor is assumed to be anisotropic so that the diagonal components are not equal to each other and the diagonal component R 12 in (9) does not vanish. While it is unlikely that an anisotropic Reynolds stress would actually develop in a homogeneous turbulent flow, this is still a useful special case in which to examine the performance of the inverse procedure before going to fully inhomogeneous, anisotropic turbulence in Subsection II B 3. The tests were performed using a rectangular domain with side lengths L x = 4 and L y = L z = 2 on a computational grid with 128, 64, and 64 points in the x-, y-, and z-directions, respectively. The computations assumed triply periodic boundary conditions, which were enforced by including one period of the SVS vortex structures in each direction, including the diagonal directions, in the velocity computation as vorticity sources. The computed Reynolds stresses were averaged over all computational points and over 10 different runs with different random vortex positions.
A listing of different prescribed Reynolds stress values used for the validation tests for homogeneous turbulence is given in Table II In each case, we plot the ratio R ij /q for each non-zero Reynolds stress component, with the predicted components on the y-axis and the prescribed components on the x-axis. The turbulent kinetic energy q was computed separately from the prescribed and predicted diagonal components of the Reynolds stress as The predicted Reynolds stresses shown in Fig. 1 are in good agreement with the prescribed values for all cases examined,  Table II. demonstrating success of the inversion procedure described in Sec. II A for homogeneous turbulence.

Validation for inhomogeneous turbulence
In actuality, anisotropic turbulence tends to form under non-homogeneous turbulent flow conditions. In this section, we examine the performance of the SVS concept and of the inversion procedure described in Sec. II A for two examples of inhomogeneous turbulent flows. In both cases, the mean flow is two-dimensional so that the Reynolds stress has the form (9), and the flow is assumed to be periodic only in the xand y-directions. The tests were performed using a rectangular domain with side lengths L x = 4 and L y = L z = 2 on a computational grid with 128, 64, and 64 points in the x-, y-, and z-directions, respectively. The predicted Reynolds stress values were averaged over the x-y plane for each z value, as well as over 20 repeated runs with different vortex positions.
The first test was for a case with isotropic prescribed Reynolds stress (R 11 = R 22 = R 33 , R 12 = 0) which varies as a top-hat distribution in y, as shown by the solid black line in Fig. 2(a). Comparison of the predicted Reynolds stresses with the prescribed distribution illustrates the nonlocal characteristics of the SVS method. The normal components of the predicted Reynolds stresses, plotted using the three color lines in Fig. 2(a), appear similar to a diffused (or filtered) form of the original profile. All three normal components are close to equal for the predicted Reynolds stress, and the predicted off-normal (shear) component (R 12 ) is close to zero.
The second test was for a case similar to an idealized turbulent planar jet, with R 11 = R 22 = R 33 and R 12 0, as shown by the solid and dashed black curves in Fig. 2(b) for the normal and shear Reynolds stresses, respectively. For simplicity, the normal stresses were prescribed as a quadratic function of y and the R 12 component was prescribed as one period of a sine wave. The predicted Reynolds stresses are in very good agreement to the prescribed values, although there is observed to be a slight flattening of the peak normal Reynolds stresses in the predicted values as compared to the prescribed values.

III. COMPUTATIONAL METHODS USED FOR VALIDATION TEST
Validation tests of the SVS method for the transport of interacting particles in anisotropic, inhomogeneous turbulence are reported in Sec. IV for turbulent planar jet flows. The current section briefly describes the computational methods used for direct numerical simulation of the fluid flow and for simulating interacting particle transport in the validation computations.

A. Computational method for direct numerical simulations
Direct numerical simulations (DNS) of turbulent planar jet flows were used to validate the SVS predictions. The DNS computations were performed using a fractional-step method, 32-34 with time advancement performed using a thirdorder Runge-Kutta method for convective terms and the 2nd order Crank-Nicolson method for viscous terms. Algorithms for all spatial derivatives except the convective terms are approximated using second-order centered finite differences (three point stencil) on a non-staggered grid. The discretized equations for the kth Runge-Kutta step are given bỹ where u and p are the fluid velocity and pressure and the coefficients α k , γ k , and ζ k are given by Rai and Moin. 32 Continuity is enforced by a projection method leading to Eq. (14c) for the pseudo-pressure, denoted by φ. In the multigrid solution of this equation, the five-point stencil produced by successive application of the gradient operation followed by the divergence operation was employed rather than a numerical approximation to the Laplacian. The Crank-Nicolson method was used to solve the Helmholtz problem, given in (14b). A tenthorder approximation was used for the convective terms, requiring an 11-point stencil. To control non-linear instabilities, at the end of each time step, the velocity components were filtered using a tenth-order filter (again using an 11-point stencil). 35,36 After filtering to obtain u filtered , the velocity u was replaced by (1 − q)u + qu filtered , with q = 0.05. The mean flow was initialized in the x-direction with cross-directional variation in the z-direction. A very weak initial turbulence was introduced using a synthetic turbulence generator, similar to the work of Smirnov et al., 37 with initial turbulent kinetic energy of 10 5 . The turbulent flow was assumed to be periodic in the xand ydirections, and a symmetry boundary condition was imposed in the z-direction. A layer of five ghost points in each direction surrounded the computational domain so that no adjustment of the differentiation schemes was needed near the domain boundaries.

B. Discrete element method for particle transport
Particle transport and collisions were computed in both the DNS and SVS computations using a soft-sphere discreteelement method (DEM) for a set of N p colliding non-adhesive particles of finite diameter d and mass m. The computations evolve the particle velocity v and rotation rate Ω by the solution of the momentum and angular momentum equations for individual particles, given by where m = π ρ p d 3 /6 and I = (1/10) md 2 are the particle mass and moment of inertia. The momentum and angular momentum equations include fluid-induced forces and torques on the particle (F F and M F ) and forces and torques resulting from particle collision (F A and M A ). The computations employ a multiple-time step algorithm to accurately resolve numerical stiffness problems introduced by the different time scales associated with the fluid flow, particle transport, and particle collisions. The time steps, including the fluid time step ∆t = O(L/U), the particle time step ∆t p = O(d/U), and the collision time step ∆t c = O(d(ρ 2 p /E 2 p U) 1/5 ), satisfy ∆t > ∆t p > ∆t c , where L and U are the characteristic length and velocity scales of the fluid flow. Here, ρ p and E p are the particle density and elastic modulus.
The fluid velocity u was interpolated from the Cartesian grid onto the particle locations with cubic accuracy using the M4 method of Monaghan. 38 The dominant fluid-induced force is the drag force, given by the Stokes drag law for low particle Reynolds numbers as Particle rotation relative to the fluid gives rise to a torque acting on the particles, where ω is the local fluid vorticity vector. Additional fluidinduced forces included in the computation include both the Saffman and Magnus lift forces, 39,40 added mass force, and pressure gradient force, as discussed by Maxey and Riley. 41 The collision forces and torques include the normal Hertzian elastic force F ne n, the normal dissipative force F nd n, the force and torque resulting from resistance to sliding [F s t S and aF s (n × t S ), respectively], and a torque M t n associated with resistance to twisting, where a denotes the particle radius. The unit normal vector n is defined by n = (x j − x i )/|x j − x i |, where x i and x j are the centroids of particles i and j, and the unit vector t S indicates the direction of relative motion of the particle surfaces at the contact point projected onto the contact plane. The Hertzian expression 42 for elastic normal force of two colliding particles is where the particle overlap δ N = a i + a j − |x i − x j | is written in terms of the radii a i and a j of particles i and j. The nonlinear spring coefficient K can be expressed as K = (4/3) E √ R, where the equivalent radius R and elastic modulus E are defined by Here, E i and E j are Young's moduli and σ i and σ j are Poisson's ratios of the two particles. The normal damping force F nd is approximated by where v R = v i − v j is the relative particle velocity, v i and v j are the particle centroid velocities, and the normal damping coefficient η N can be related to the restitution coefficient e using an expression due to Tsuji et al. 43 The current computations are performed with a fixed, small value of the restitution coefficient (e = 0.10), which is consistent with the observation that particle collisions occur in this problem with small values of the Stokes number, St = ρ p d 2 U/18 µL. A spring-dashpot-slider model is used to approximate the sliding resistance. 44 In this model, the sliding force F s is first absorbed by the spring and dashpot until its magnitude reaches a critical value F crit = µ f |F n |. The friction coefficient µ f is selected to have a value of 0.3, which is in approximately the middle of the range of typical values for dry surfaces discussed by Johnson. 45 If |F s | > F crit , then the colliding particle surfaces slip relative to each other and the friction coefficient is given by the Amonton expression, For the subcritical case |F s | < F crit , the sliding resistance due to the spring and dashpot for particle i is given by where the slip velocity v S (t) is defined by and the slip direction is t S = v S /|v S |. The time integral in the first term in (22) gives the tangential elastic displacement of the material before slipping occurs, where t 0 is the time of initial particle impact. The expression for the tangential stiffness coefficient k T derived by Mindlin 46 can be expressed as where G −1 ≡ 2−σ i G i + 2−σ j G j is the equivalent shear modulus and G i = E i /2(1 + σ i ) and G j = E j /2(1 + σ j ) are the shear moduli of the two particles. We follow Tsuji et al. 43 in assuming that the tangential dissipation coefficient is of the same order as the normal viscous damping coefficient, and thus set η T = η N .
Twisting occurs when the two colliding particles have different rotation rates in the direction n. The relative twisting rate Ω T is defined by The twisting resistance force is given by where the time integral represents the angular displacement prior to torsional sliding. Expressions for the torsional stiffness and viscous friction coefficient are similarly given by Ref. 47, The particles begin to spin relative to each other when the torque exceeds a critical value, given by When |M t | > M t,crit , the torsional resistance is given by

A. Direct numerical simulation
Direct numerical simulations were conducted for a particulate turbulent planar jet flow with one-way coupling between the fluid and particles, and the results were compared to SVS simulations of the same problem. The computational domain was discretized using a Cartesian grid over the intervals −2 ≤ x/H ≤ 2, −3 ≤ y/H ≤ 3 and −2 ≤ z/H ≤ 2, where 2H is the initial jet width. The use of a uniform grid with (N x , N y , N z ) = (129, 193, 129) points for DNS led to grid increments that were nearly the same in all directions. The initial jet Reynolds number is given by Re H = U 0 H/ν = 3200, where U 0 is the nominal jet velocity and ν is the kinematic viscosity. The initial mean velocity profile U(z) was chosen to be of the hyperbolic tangent form where θ 0 is the initial momentum thickness and the centerline of the jet corresponds to y = 0. For the current computations, we select H/θ 0 = 35, for which value da Silva and Pereira 48 showed that the most unstable Kelvin-Helmholtz instability wavelength is λ KH = 0.87H, which is less than onequarter the grid domain length in the streamwise direction. The DNS simulations were performed using a fixed time step of ∆t = 0.005H/U 0 , which was selected to yield a Courant-Friedrichs-Lewy (CFL) number less than 0.1. In the following, length, velocity, and time are nondimensionalized by H, U 0 , and H/U 0 , respectively. Results for time variation of the turbulent kinetic energy q, the dissipation rate per unit mass ε, and enstrophy per unit volume Ω for the jet flow are shown in Fig. 3, defined by where D ij are the components of the rate of deformation tensor, u and ω are the velocity and vorticity vectors, respectively, and V ave is the averaging volume. Since we want these measures to be independent of the size of the computational domain, we performed the averaging only over the region −1 ≤ y ≤ 1 initially occupied by the jet. The turbulent kinetic energy initially increases as the turbulence develops in the jet, up to a time of t 10.5, at which the peak value of q is observed. It then gradually decreases as the turbulence within the jet decays. The time variation of the dissipation rate and enstrophy also exhibits an increase at the beginning of the computation, a peak, and then a gradual decrease, although the peak values for enstrophy and dissipation rate occur a little later than for kinetic energy (close to t ≈ 12). Contours of the velocity magnitude at a series of times during the jet development are shown for both the DNS computation and the comparison SVS computation over a series of times in Fig. 4. These contours illustrate the development of instability and turbulence at the beginning of the run (t ≤ 10), followed by decay of both the mean jet velocity and the turbulence within the jet at later times (t ≥ 15). The turbulence decay is accompanied by outward spreading of the turbulent velocity field and decrease in velocity magnitude values within the central region of the jet. The jet decay is often characterized in the similarity theory by two time-varying parameters-the centerline velocity U c and the distance δ 1/2 from the centerline at which the mean velocity equals one-half the centerline value. The former of these parameters characterizes the jet strength, and the latter characterizes the jet width. It is recalled that in their experiments with a spatially varying planar jet, Gutmark and Wygnanski 49 observed that δ 1/2 and 1/U 2 c both vary approximately linearly with distance. This observation suggests that by replacing the downstream with time might be observed for the current problem of a temporally varying jet. Figure 5 plots time variation of both δ 1/2 and 1/U 2 c , exhibiting nearly linear variation in time in both cases.
A comparison of the mean velocity and normal Reynolds stresses from our DNS computations with results from previous experimental and numerical studies is given in Fig. 6. The values are non-dimensionalized using U c (t) and δ 1/2 (t) to write them in similarity form, and we have confirmed that the results are nearly independent of time during the similarity regime of the computation (10 ≤ t ≤ 20). The mean velocity curve from our DNS results is very close to the mean velocity in the comparison studies. The normal Reynolds stress results are also reasonably close to the values in the comparison studies, although the Reynolds stresses exhibit more scatter among the different studies than do the mean flow results.

B. Comparison of DNS flow field to SVS results
The grid used for the SVS computations had (N x , N y , N z ) = (128, 128, 128) points. It is a requirement of the accelerated method used for the SVS method that the number of points on each side be a multiple of two. The SVS simulations were conducted using the DNS Reynolds stress results over the interval 10 ≤ t ≤ 20 for which the similarity solution was found to be valid in the DNS results. Before this time period, the DNS results show that the turbulence is still developing, and after this time period, the turbulence exhibits rapid dissipation. The Reynolds stress predictions from DNS were written in dimensionless similarity form (as shown in Fig. 6) and averaged over the computational time period 10 ≤ t ≤ 20, in order to smooth out temporal fluctuations. These averaged Reynolds stresses in similarity form were then read into the SVS simulations, along with the DNS predictions for U c (t) and δ 1/2 (t) shown in Fig. 5, and used to generate time-varying prescribed Reynolds stress profiles for use during the SVS computation. The SVS computation was initialized with a prescribed number N V of vortices positioned randomly in the SVS domain. The vortex strength and orientation were set using the prescribed Reynolds stress field at t = 10, obtained from the DNS results as described earlier, using the inversion method described in Sec. II A. A plot showing the initial strength distribution and initial orientation of the SVS vortices is given in Fig. 7. While the vortices were located throughout the computational domain, the vortices with significant strength were located primarily within the interval −1 ≤ y ≤ 1. All initial vortices were randomly assigned an initial "age," which advanced with time during the computation. When a vortex age exceeded the prescribed vortex lifespan T V , the vortex was removed and a new vortex was introduced at a random location within the computational domain. The strength and orientation of the new vortex were again set using the procedure described in Sec. II A using the prescribed Reynolds stress field for the time that the vortex is introduced. Consequently, as the turbulence decays in time during the SVS computation, the strength of the newly initiated SVS vortices generally decreases at a given position in the flow field. A series of SVS computations with different values of N V and T V were performed, as listed in Table III. The "standard" SVS computation (case S) was selected as one with N V = 1024 and T V equal to the integral time scale T 0 at t = 10.
A comparison of the time variation of the velocity magnitude contours for the SVS generated flow field at t = 10, 15, and 20 is given in Fig. 4 immediately below the DNS plot at the same time (and using the same color scale). We do not expect exact agreement since the SVS vortex structures are randomly distributed in space, but it is noted that the velocity magnitudes and general tendencies of the SVS generated flow field are similar to the DNS flow. In both cases, the simulated jet turbulence gradually spreads in the y-direction and decays over this time interval. As would be expected from the uniform vortex size used in the SVS formulation, we observe that the DNS flow field results in Fig. 4 exhibit more small-scale structures than do the SVS flow fields. A plot showing the time variation of the jet centerline velocity U c and the jet width measure δ 1/2 is given in Fig. 8. The value of U c decreases during the time interval and the value of δ 1/2 increases, as expected for decaying turbulence. The SVS predictions for U c and δ 1/2 are observed to be significantly noisier than the DNS predictions. This noise in the SVS predictions is associated with the "death" of some vortices and the "birth" of new vortices at random positions in the flow field. The SVS predictions for U c fluctuate closely about the DNS results. The SVS predictions for δ 1/2 are also close to the DNS predictions in the beginning part of the computation (t ≤ 14), but by the end of the computation, the predicted jet width measure for SVS is about 10% lower than that for DNS. Similar fluctuations in the SVS predictions are shown in Fig. 9(a), in which we compare the time variation of the turbulent kinetic energy for the DNS and SVS computations. The SVS result is again observed to fluctuate around the smoother DNS prediction, with a root-mean-square value that decreases when the value of the vortex lifetime T V is reduced. The power spectrum is plotted in Fig. 9(b) at time t = 15 for both the DNS and SVS computations. Both computations exhibit a k −5/3 Kolmogorov spectrum in the inertial range, with DNS and SVS spectra in close agreement. At a high wavenumber, the SVS spectrum reduces much faster than the DNS spectrum as a consequence that SVS contains only vortices with the length and velocity scaled to the integral scale eddies.
A comparison of the time-averaged Reynolds stresses, nondimensionalized using the similarity variables, is given for  Table III. DNS and SVS in Fig. 10. The DNS values of U c and δ 1/2 are used to write the Reynolds stresses and lateral distance in similarity form for both computations. The three normal Reynolds stress values are very close for the DNS and SVS predictions. The SVS prediction for the dimensionless Reynolds shear stress R 12 /U 2 c exhibits lower peak values than the DNS predictions but otherwise has the same form.

C. Comparison of DNS particle transport to SVS results
An initial DNS flow computation was conducted out to a time of t = 10 with no particles in order to allow the turbulence to develop and to achieve a self-similar state. The DNS computation was then restarted with particles present and continued out to a time t = 20. A total of N p = 32 000 particles of diameter d = 0.04 and density ratio χ = ρ f /ρ p = 1 were used. The particles were initially placed randomly within the region −1 ≤ y ≤ 1 covering the jet. The particle Stokes number based on the jet width scaling, St H , is given by Particle initial positions were identical for both the DNS and SVS simulations. The particle concentration profile in y was computed by dividing the flow field into bins and then adding the volume of particles contained in each bin. For particles that straddle the boundary between bins, the particle volume is divided along the bin boundary and only that portion of the volume lying in each bin is included in the sum. The concentration field for SVS and DNS is identical at the initial time t = 10 and has a top-hat form as shown in Fig. 11(a). During the time period of the flow computation 10 ≤ t ≤ 20, the concentration field spreads outward into the lateral regions around the jet due to forcing by the jet turbulence. The resulting concentration field for both DNS and SVS computations at time t = 20 is plotted in Fig. 11(b), exhibiting excellent agreement between the two methods. This comparison demonstrates that the SVS method accurately simulates dispersion of the particle field. Another way to examine particle dispersion is to calculate the root-mean-square particle position y rms , defined by where y n,par denotes the y-position of particle n. A comparison of y rms as a function of time for DNS and for a variety of SVS computations with different parameter values is plotted in Fig. 12. Figure 12(a) shows the effect of number of vortices N V on the lateral particle dispersion in cases with T V = T 0 .
As the number of vortices decreases, the strength of each vortex is increased so as to hold the turbulent kinetic energy fixed. As can be seen, cases with smaller number of vortices (e.g., N V = 256) exhibit slower lateral dispersion, resulting in lower values of y rms at the given time than the DNS predictions. At higher number of vortices, the predictions of the various SVS computations appear to converge to a value of y rms that is close to the DNS prediction up to a time of about t = 17, after which the SVS predictions are somewhat less than that for DNS. Figure 12(b) shows the effect of vortex lifetime on lateral particle dispersion. Increase in vortex lifetime is found to increase the rate of particle dispersion from the center of the jet, up to a lifetime value of about T V = 1.5T 0 , above which the particle dispersion rate remains close to the DNS prediction. This increase in dispersion rate occurs because longer residence of strong vortices near the jet center allows them longer time to repel particles via centrifugal force. We also note that the turbulent kinetic energy in the SVS computation increases (above the DNS prediction) as the vortex lifetime is increased significantly above the integral time scale T 0 , which also increases the lateral dispersion rate. The total number of particle collisions is plotted as a function of time in Fig. 13 for DNS and for a variety of SVS computations with different values of N V and T V . Figure 13(a) shows that the number of collisions in SVS computations is lower than that for DNS for small numbers of vortices but that the collision number increases to close to the DNS results as FIG. 11. Particle positions (a) at the start of the particle runs at t = 10 and (b) at the end of the run at t = 20 for DNS (red) and SVS case S (blue). (The particle positions at t = 10 are the same for DNS and SVS.) the number of vortices increases. The variation of vortex lifetime is seen in Fig. 13(b) to have little effect on the number of particle collisions, which we believe to be a consequence of two opposing influences. As discussed previously, increasing the vortex lifetime tends to disperse the particles more rapidly in the y-direction, consequently decreasing particle concentration and leading to lower numbers of collisions. On the other hand, increasing the vortex lifetime also introduces a lag that increases the turbulent kinetic energy slightly in a decaying turbulent flow, resulting in an increase in number of particle collisions. These two phenomena counteract each other so that little change in collision number with vortex lifetime is observed in Fig. 13(b).
The tendency of particles to cluster can be characterized by the radial distribution function (RDF), g(r), which is 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 obtained by computing the average number of neighboring particles whose centroids are located within a radial distance r from a given particle centroid. In order to smooth the RDF values, we have averaged the predicted RDF for both DNS and SVS over the time interval 14 ≤ t ≤ 16, which was selected because this time interval is in the middle of the computational interval (10 ≤ t ≤ 20). It is sufficiently small that the turbulence kinetic energy does not change by a large amount, and yet it is also sufficiently large that noticeable smoothing of the data is observed. The radial distribution function is plotted in Fig. 14 for both DNS and SVS computations and found to compare well. The RDF peak in the SVS computations is a little higher than the DNS result, which might be a consequence of the observation that DNS was observed to disperse particles a little more rapidly in the lateral y-direction, and so the resulting concentration is slightly lower, but the effect is small.

V. CONCLUSIONS
The paper presents a novel inverse method by which the orientation and strength of a set of finite-length vortices can be set to reproduce a prescribed anisotropic Reynolds stress field. This inverse method was incorporated into the stochastic vortex structure (SVS) algorithm to generate a time-varying synthetic turbulence field for transport of interacting particles in anisotropic, non-homogeneous turbulent flows. The proposed SVS method is well suited for simulation of interacting particles since the statistics of the generated synthetic turbulence is both structurally and temporally consistent with the original turbulence and it can be computed rapidly with use of the fast multipole accelerated method. 28 It has been previously demonstrated 27,28 for homogeneous, isotropic turbulence that the SVS method accurately reproduces the turbulence energy spectrum, the probability density function of the acceleration, velocity and vorticity fields, the collision rate of advected particles, and a variety of agglomeration measures (fractal dimension, size distribution, etc.) for adhesive particles. The current paper extends the SVS approach to make it a viable method for arbitrary turbulent flows and not only for homogeneous turbulence.
The effectiveness of the proposed inverse method was demonstrated in a series of computational experiments. We first examined the accuracy of the inverse method for an anisotropic, but homogenous, turbulent field with different prescribed values of the Reynolds stresses. Next, we examined the performance of the inversion procedure for setting the initial vortex orientation and strength in two different nonhomogeneous turbulent shear flows. Prescribed and predicted Reynolds stresses were compared for the above cases and showed good agreement. Finally, the SVS predictions for flow and particle transport in a planar turbulent jet flow were compared with direct numerical simulation (DNS) results. The SVS computations used the Reynolds stress profiles computed from DNS together with our inverse procedure to specify the initial orientation and strength of the stochastic vortices, both at the start of the computation and when new vortices were introduced during the computation. The Reynolds stress profiles of both DNS and SVS computations were normalized in similarity form and averaged over the duration of the SVS computation and found to compare well. Measures of particle dispersion, clustering, and collision during the SVS and DNS computations were also found to be in good agreement. The effect on the SVS predictions of variation of the number and lifetime of vortices was also investigated, as these are two important numerical parameters that must be specified in the SVS computations. Computations with small numbers of vortices yield too low collision rate and weak dispersion, but the results approach the DNS predictions as the number of vortices is increased. The particle dispersion predictions were poor when the vortex lifetime was significantly below the turbulence integral time scale, but values near the integral time scale up to about twice the integral scale yielded acceptable