Algorithms and their Applications
One of the perennial issues with numerical codes, especially those geared towards a large amount of "number crunching" is that they require unfavourable amount of wall-clock time. This section presents a brief discussion of the physical problem, the algorithms used and some sample results. Particularly discussed are,
-
k-means clustering - Group particles together in brownout calculations to reduce the cost
-
Depth first search - Detect regions of vorticity in highly congested flow fields
-
kd-tree - Search algorithm to find the wall distance in CFD meshes
-
Fast Multipole Method and oct-trees - Efficient method to solve N2 interaction problems (vortex/vortex interaction)
-
Panel methods - Predict the effect of blade erosion
-
GMRES - Quick method to invert sparse matrices for typical implicit solvers in CFD (currently in development...)
k-means Clustering (Track particles in brownout clouds)
The motivation stems from the requirement to track a prohitively large number of dust particles in brownout calculations. The k-means clustering method is based on the principle that certain sets of individual particles can be decomposed into smaller groups of clusters, and that the resulting equations of motion are solved only for the clusters. The displacements of the cluster centers are then applied to all of the individual particles comprising the cluster. There are two aspects to this method. 1. Selecting the candidate particle groups to form a cluster; 2. The effects produced by clustering, declustering, and reclustering of the particle groups.
The general idea of the k-means clustering approach is depicted in Fig. 1. The k-means clustering algorithm is a category of cluster analysis that aims to partition n observations (i.e., particles) into k clusters in which each particle belongs to the cluster with the nearest mean. The most common algorithm uses an iterative refinement technique.
More on the application of this algorithm to study brownout can be found here.
Fig. 1: Schematic showing the k-means method of clustering and declustering.
Depth First Search (Detect vortices in congested flow fields)
One of the quantities of interest in aerodynamic flows is that of vorticity and circulation. Particle image velocimetry (PIV) techniques typically capture snapshots of the velocity in the field from which, the vorticity is derived and the circulation is obtained using the box circulation method. The problem with the box circulation method is that the vortex or structure of interest must be isolated so that the box drawn does not cross extraneous flow structures. Finding a "box" robustly can become increasingly challenging for complicated flows, such as the trailed wake sheet off of a hovering rotor or eddies in the boundary layer..
A depth first search algorithm was used the compute the circulation from the vorticity field based on Stokes theorem (vorticity was obtained from the velocity field). The advantage of this algorithm is primarily two-fold: 1. No requirement is placed on the shape/size of the flow structure, and 2. The search is of O(N), where N is the size of the image. The search algorithm sweeps through the image searching for pixels (or grid cells) that belong to a vortex blob and eventually finds all such vortex blobs. Figure 2 shows the flow field near the blade tip at a wake-age of 3 degrees, with both positive and negative vorticity identified. Figure 3 shows the circulation at the tip and that compared against a strip theory that is typically used in rotorcraft flows (BEMT - Blade Element Momentum Theory).
Fig. 3: Extracted circulation profile using the depth-first search.
Fig. 2: PIV data of vorticity field beneath a hovering rotor.
kD-tree Search (Find wall distance for CFD turbulence models)
Finding the wall-distance is essential to most turbulence modelling in CFD applications. This distance is defined as the distance of a volume cell (field point) to the nearest solid boundary surface (wall point). If we assume for simplicity that there are N field points and N wall points, the cost of find the wall points using a naive search is O(N^2).
A kD-tree was used to greatly simplify the cost of finding the nearest wall-distance. A kD-tree is essentailly a binary tree with each entry being of possible k dimensions. A graphical representation of the wall points and a representative field point is shown in Fig. 4. In the current implementation, the elements are sorted using a heapsort and the median is inserted into the tree. This procedure results is a well-balanced tree, required for efficient searching. As shown in the results in Fig. 5 the kD-search allows for O(N log^2 N) search and typically only 60-70 points are visited for N = 1,000,000 with the time dropping from 2.7 hours to 2.3 seconds.
Fig. 5: Comparison between the naive search and kD-tree search.
Fig. 4: Schematic representation of the kD-tree.
3D Multi-Level Fast Multipole Method (N-body interactions)
Regarded as one of the top ten algorithms of the previous century, the Fast Multipole Method (FMM) is an extremely efficient method to solve for matrix-vector product wherein the entries of the matrices (kernel function) have a specific form. Typical applications are those that require the solution to N-body problems; such as intereactions between heavenly bodies, electric charges, and vortices. If we think about gravitational forces, the effect of a cluster of planets or a galaxy far away on a planet/star can be simplified by assumed one mass for the galaxy, i.e., it is not required that the effect of the individual constituents of the galaxy be present. This idea of splitting the domain into near-field and far-field forms the basic idea behind FMM.
An oct-tree spatial partitioning data structure was used to split the domain in 3D, which allows for log N traversal up and down the tree. The S|S, S|R, and R|R translations were built from the spherical harmonics and the direct evaluations were computed using the Biot-Savart kernel in the present problem of interacting vortices (Fig. 6). More information on FMM can be found in this seminal paper by Greengard and Rokhlin.
Figure 7 shows the current state of the code and shows the CPU time as ut varies with N. While the expected O(N log N) efficiency was obtained, the breakeven point is very large (~N = 20,000). The reason behind the large breakeven point is the lack of precomputed translation stencils and the need for a "point and shoot" algorithm for translation operators.
This algorithm is still under development.
Fig. 7: Preliminary results obtained using MLFMM.
(The expected O(N log N) exists but breakeven can be reduced)
Fig. 6: Elements of FMM; Spherical harmonics, oct-tree data structure and the order of translations from source to target.
Panel Methods (Study and predict blade erosion)
Blade erosion occurs because of the interaction of the rotor blades with the suspended particulates causing wear and tear of the leading edge of the blades. In a brownout environment, this erosion process is greatly accelerated. To prevent damage from dust/sand erosion, metallic abrasion strips, typically composed of stainless steel or titanium, are bonded to the blade leading edge to serve as a hard surface that absorbs the kinetic energy of the sand particle. These impacting sand particles, however, slowly erode away the metallic material. Current efforts include preventive measures such as coatings and abrasion strips to delay the erosive effects of blade erosion.
A 2D panel method was developed (a 3D panel method was also developed and was used here) and the particles were tracked based on the velocity induced on it from the panels and their collisions with the airfoil surface was subsequently modelled. Interesting physics were observed and estimates of the energy transferred to the airfoil were also obtained; see Figs. 8-11.
Fig. 9: Schematic of the process of a particle colliding with the airfoil surface.
Fig. 11: Trajectories of particles of different diameters as they collide with a NACA 0012 airfoil.
Fig. 10: Kinetic energy loss per particle of different sizes upon collision with an airfoil.
Fig. 8: Validation of the panel method.