|
The capablity of vortex methods to efficiently solve problems in incompressible fluid flow is well-known. Two difficulties of modern vortex methods are numerical dissipation of vorticity and lack of results for long-term runs. In my research, I attempt to address these issues by using a persistent, connected, vortex sheet discretization method, and VIC method.
On 2006-04-27, I successfully defended my dissertation. You may download a PDF copy with the following link: mstock_dissertation.pdf (7.0MB).
In the process of conducting this research, I have assembled into one document a summary of the major aspects of vortex methods research as well as 415 literature references. The document is available for free here: vortex_methods_literature.pdf (723kB, 2007-04-18), and you can also grab the BibTeX file for it. It is a living document, and thus contains not only references and questions, but my throughts regarding various lines of research. I am making it available in the hopes that newcomers to the field might find it useful.
The majority of computational fluid dynamics calculations use an Eulerian description of the flow, in which the governing equations are written as finite difference equations relating properties at two or more fixed and adjacent points in space. The equations are solved at every time step for the new properties at each point. In contrast, a Lagrangian description of a flow involves tracking the motion of fluid parcels, each of which is allowed to carry flow properties with it. This generally makes the advection term, which commonly causes unwanted dissipation in Eulerian schemes, simple.
Rosenhead (1931) was the first to represent a flow as a collection of potential flow vortex elements. The flowfield was calculated as the superposition of the flowfields of the individual vortex singularities. One advantage of this method is that conservation of mass is always satisfied. Because vorticity can be more compactly represented than velocity, complex flows can be simulated with a smaller number of elements than are required for an Eulerian description.
The vorticity in a flow can be represented by independent particles, connected filaments, or connected sheet elements. Particles are the easiest to implement, because no connectivity information needs to be stored, and remeshing to control growth of the region of vorticity is easy (but dissipative). In a vortex particle method, each particle has an amount of vorticity and a volume assigned to it. The resulting vorticity field is simply a sum of the contributions of these individual vorticities. Three disadvantages to the vortex particle method are apparent, though. First, in 3D flows, the vortex stretching term must be taken into account in a conservative manner, which consists of an extra step and careful formulation of the correction term. Second, remeshing, which is neccessary to fill in the voids left as the particles advect, causes numerical diffusion. Finally, incorporating multiple fluid phases in the method is difficult.
A vortex filament method discretizes the vorticity into short segments of vortex lines of constant circulation. Naturally, the connectivity must be stored, but thanks to Kelvin's circulation theorem, the circulations stored with the filaments do not have to be updated to correct for the vortex stretching term. Vortex filaments, though, still cannot adequately define the boundary between two fluid phases without complex remeshing strategies.
Vortex sheet methods can use connected triangles, quadrilaterals, or other smoother representations of a surface, upon which a vector-valued velocity jump, vortex strength, or circulation is carried. These methods, too, must store connectivity information, and, like filament methods, can leverage Kelvin's circulation theorem to bypass the vortex stretching term. A vortex sheet is also well-suited to define the boundary between two fluid phases, and correction terms can be applied to the strength of the vortex sheet if the phase boundary also represents a density discontinuity. One major disadvantage, though, of both filament and sheet methods is the difficulty dealing with the exponential growth of the area of the vortex sheet (or length of vortex lines).
1.3.1 Direct summation. The direct solution method consists of integrating the Biot-Savart equation completely, which requires a separate calculation for every pair of participating vortex elements. The number of operations is obviously order N^2, making this method unattractive for any real engineering calculation. Additionally, it is believed that singularity formation in finite time may be unavoidable without regularization. Krasny (1986) introduced regularization to this method by desingularizing the kernel, effectively turning the point vorticies into finite-core vorticies.
1.3.2 Vortex-In-Cell. (Christiansen, 1973) A variant of the particle-in-cell method created for electrostatic calculations, the Vortex-In-Cell, or VIC, method uses a temporary grid to provide regularization and to speed the method. Using fast Poisson solvers, the method achieves a computational cost of order M log M + N, with M being the number of cells in the temporary grid, and N being the number of vortex elements.
1.3.3 Fast multipole method. (Greengard and Rokhlin, 1987) This recent development in multiresolution modelling uses the spherical harmonic expansion of the potential function to allow more accurate particle-cluster interaction computations. Allowing box-box interactions, the computational cost of FMM scales as order p^2 N, where p is the order of the multipole expansion, and N is the number of particles. Simulations with millions of particles have been run with this method. (Koumoutsakos and Leonard, 1995; Winckelmans and others, 1995)
1.4.1 Species interface/mixing. Using methods standard in front-tracking simulations, a three-dimensional front (defined by triangles or quads, etc.) can define a boundary between fluid phases, or of a marker fluid. Some methods include a term to account for surface tension. Additionally, in the case of an interface between miscible fluids, molecular or subgrid-scale dissipation can be modelled by allowing each surface element to carry a series of scalar moments, together describing the surface-normal profile of the scalar gradient. These scalar moments can be updated in time to reflect the effects of shear on the thickness of the mixing interface. Care must be taken not to allow the mixing interface to become too thick, as that would introduce errors in the advection term. Remeshing the single layer into a series of layers alleviates those errors.
1.4.2 Density interface/stratification. This boundary surface, which can concurrently be a vortex sheet and have scalar gradients defined as in section 2.4.1, can also carry a density gradient. In the case of small density gradients, characterized by an Atwood number close to zero, simplifications of the equations allow trivial calculation of the vortiticy generated on the boundary surface.
1.4.3 Large Eddy Simulation. The regularization imposed by some vortex methods can be seen as a sort of subgrid-scale model. For example, element or filament surgery has a subgrid dissipative effect. In order to create an accurate subgrid- scale model, though, one would not only have to quantify or control the effects of regularization, but also add an appropriate term to account for the transfer of energy (or enstrophy) to and from the unresolved scales not accounted for by regularization. Viscous vortex schemes have been created to more accurately treat the diffusion term in the Navier-Stokes equations, most notably the Particle-Strength-Exchange scheme (Mas-Gallic, 1987). Cottet (1996) used a PSE formulation to model the subgrid-scale energy transfer in a Lagrangian vortex method.
The vortex sheet is a three-dimensional triangulated mesh, with each element pointing to three nodes, and each node being referenced by all adjacent triangular elements. Each element stores three circulations, corresponding to its three edges. These circulations uniquely describe a vector-valued vortex sheet strength for the element. The advantage of this method is that the elements' edge circulations need not be corrected for vortex stretching. The circulations are only modified due to baroclinic generation of vorticity, or due to a dissipation model.
The code begins the velocity calculation stage as a standard VIC method would--by placing all of the elements' vorticities onto a temporary grid. This can be done with trilinear interpolation (volume-weighting), rectangular M4', rectangular Peskin, or conserving spherical Peskin functions. The negative of the curl of the vorticity field (with appropriate boundary treatment) becomes the right-hand side of a Poisson equation for velocity. This method is different from previous methods that solve the Poisson equation for a three-dimensional streamfunction. The temporary vorticity field created in this step is used only to calculate the velocities, and is discarded afterward. Currently, the FISHPAK solver is used to solve the Poisson equation for the velocity field.
The vortex sheet strength is modified at each time step to account for baroclinic generation of vorticity. This change is proportional to the area of the element and the Atwood number, and is in the direction of the cross-product of the element's surface normal and a term consisting of gravity and local accelerations. The vector-valued result is converted to a trio of edge circulations and then added to the element.
To maintain the order of the method despite the continual stretching and folding of the interface, several remeshing steps are undertaken. Maintenance of element (triangle) size is accomplished by finding any element edge longer than a threshhold and splitting the triangles sharing that edge. The new node is placed at the midpoint of the long edge, and a proper amount of vortex sheet strength is transferred from the old, large triangular elements to the newly created ones.
Splitting can only limit element edge lengths, and, thus, element areas. In order to reduce the number of elements while maintaining mesh resolution, close nodes within the same sheet are detected and merged together. Again, the connecting elements' vortex sheet strengths are modified to ensure that the method is conservative.
Despite this aggressive remeshsing, long runs can still create large amounts of surface. Methods to merge overlapping sheets are currently under investigation.
Each element has, assigned to it, a scalar jump. In the scalar reconstruction step, the scalar jumps from each element are converted to vector-valued scalar gradients and placed on the grid using the same routines used for vorticity. In this case, the divergence of this gradient field is the right-hand side of a Poisson equation for the scalar fraction field. This calculation is only performed when necessary for output, and the scalar field values are not stored or reused.
Some past simulation results are shown below, ordered from oldest (top) to most recent (bottom):
|
|
|
|
|
|
|
vort3d is a vortex method flow solver written in ANSI C as part of my PhD research with Prof. Werner Dahm at the University of Michigan and Prof. Gretar Tryggvason at WPI. A commercial product using similar techniques is NGB Technologies' LIM3D.
vort3d currently includes the following features:
Current research areas include the following:
It is clear that a vortex sheet method is capable of capturing the salient features of a variety of incompressible flows, especially multiphase flows with small density differences. Future research will push the applicability of these methods to multiphase flows with large density variations, long duration runs, combustion, and Large Eddy Simulation.
J. Christiansen, Numerical simulation of hydrodynamics by the method of point vorticies, J. Comput. Phys. 13 (1973) 363-379.
G.-H. Cottet, Artificial viscosity models for vortex and particle methods, J. Comput. Phys. 127 (1996) 299-308.
L. Greengard and V. Rokhlin, A fast algorithm for particle simulation, J. Comput. Phys. 73 (1987) 325-348.
P. Koumoutsakos and A. Leonard, High resolution simulations of the flow around an impulsively-started cylinder using vortex methods, J. Fluid Mech. 328 (1995) 1-38.
R. Krasny, Desingularization of periodic vortex sheet roll-up, J. Comput. Phys. 65 (1986) 292-313.
S. Mas-Gallic, Contribution a l'analyse numerique des methodes particulaires, These d'Etat, Universite Paris VI (1987).
L. Rosenhead, The formation of vorticies from a surface of discontinuity, Proc. Roy. Soc. London Ser. A 134 (1931) 170.
M. Stock, Summary of Vortex Methods Literature (2006)
G.S. Winckelmans and others, Application of fast parallel and sequential tree codes to computing three-dimensional flows with the vortex element and boundary elements method (1995), in Vortex flows and related numerical methods II, ESAIM Proc. (1996) 225-240.