WO2001001237A1 - Apparatus for and method of simulating turbulence - Google Patents

Abstract

In accordance with a preferred embodiment of the invention, a novel apparatus for and method of simulating physical processes such as fluid flow is provided. Fluid flow near a boundary or wall of an object is represented by a collection of vortex sheet layers. The layers are composed of a grid or mesh of one or more geometrically shaped space filling elements. In the preferred embodiment, the space filling elements take on a triangular shape. An Eulerian approach is employed for the vortex sheets, where a finite-volume scheme is used on the prismatic grid formed by the vortex sheet layers. A Lagrangian approach is employed for the vortical elements (e.g., vortex tubes or filaments) found in the remainder of the flow domain. To reduce the computational time, a hairpin removal scheme is employed to reduce the number of vortex filaments, and a Fast Multipole Method (FMM), preferably implemented using parallel processing techniques, reduces the computation of the velocity field.

Description

APPARATUS FOR AND METHOD OF SIMULATING TURBULENCE

Government Interest: This invention was made with government support

under Grant No. DE-FG02-97ER82413 awarded by the Department of

Energy. The government has certain rights in the invention.

BACKGROUND OF THE INVENTION

For over a century there has been the desire and need to study and predict

physical processes such as fluid flow. Understanding the flow past an object

(e.g., automobile, airplane wing, etc.) requires an understanding of the

turbulence introduced as a result of the relative motion of the fluid and the

object. There have been many attempts in the past to gain a better

understanding of fluid flow and turbulence.

U.S. Patent No. 3,787,874, for example, disclosed a method of making

boundary layer flow conditions visible by a reactive layer of a colored chemical

indicator to the surface of an object exposed to fluid flow. U.S. Patent No.

4,727,751 disclosed a mechanical sensor for determining cross flow vorticity

characteristics. The sensors contained hot-film sensor elements which operate as

a constant temperature anemometer circuit to detect heat transfer rate changes.

These known approaches require the production of physical models and

physical testing. Such approaches are often time-consuming and expensive to

perform in conjunction with a design process as any changes suggested by the tests would require modification of the model and re -testing with the modified

model.

Computer modeling techniques that permit the simulation and analysis of

fluid flow using computer modeling is much more attractive because it eliminates

the need for physical models and repetitive testing. Some computer simulations

attempted to model the fluid flow by covering the entire fluid flow domain with

a large mesh or grid. Such grid-based simulations were inefficient and lacking in

accuracy to the real physics in fluid flows, particularly those with high Reynolds

number turbulence.

With the breakthrough work of Alexandre J. Chorin, University of

California in Berkeley, however, a tremendous advantage in understanding and

predicting fluid flow was attained using Chorin's vortex method. According to

the vortex method, the vorticity of fluid was used to better understand the basis

for fluid flow. In the vortex method a flow field is represented by a collection of

convecting and diffusing vortex elements. This method has negligible if any

numerical diffusion so that high Reynolds number turbulence effects can be

represented with high numerical accuracy and fidelity. The method is grid-free

so that it can be easily applied to engineering flows in complex geometries.

Moreover, the method is naturally adaptive since the computational elements

occupy only the relatively small support of the flow field where the vorticity is

significant.

%. The direct implementation of vortex methods, however, presents at least

two problems that have limited their appeal. First, vortex methods

indiscriminately track the dynamics of individual vortical structures as they stretch

and fold, as a result, the number of vortices can grow to impractical levels.

Second, the computational time associated with the nominally O (N_v ²)

operations implicit in the Bio Savart law evaluation of velocities due to N_v vortex

elements, an essential aspect of their formulation and a source of their unique

advantages, can grow beyond reasonable levels for engineering use. An

additional concern, which pertains to the traditional random vortex method, is

the noise introduced by a random walk diffusion model: unless a very fine

coverage of vortex elements is employed, this can easily exceed the levels found

in real turbulent flows.

The options for taking account the essential process of vortex stretching

over a field of computational elements is a strong function of the particular types

of vortices used in the numerical algorithm. The commonly employed vortex

blobs, i.e., independently evolving three-dimensional regions of high vorticity

concentration, have many desirable properties in the computation of laminar flow

solutions, yet are a source of lingering problems when turbulent conditions are

encountered. In particular, they tend to be associated with numerical instabilities

arising, apparently, when the non-linear coupling between blob strength and

stretching causes the former to grow without bound. A common tendency has been to " switch off" the vortex stretching effect to prevent instability, even

though this eliminates one of the most fundamental processes of turbulent flow.

New vorticity most often appears in viscous flows at solid boundaries as a

natural consequence of the no slip condition, by which fluid in contact with a

solid body must move with the same velocity as the solid. Under turbulent

conditions, the regions where new vortices appear tend to be relatively thin

boundary layers with very high levels of vorticity. The challenge for a vortex

method is to efficiently represent such large thin regions with vortex elements

while satisfying the no slip condition and generating the proper amount of new

vorticity. Typical vortex elements such as blobs and tubes are poorly constructed

to represent thin, flat regions, and vortex methods which attempt to use such

elements to satisfy no slip tend to be inefficient: large numbers of overlapping

elements are needed to counter the inherent three-dimensionality of the velocity

field they produce.

SUMMARY OF THE INVENTION

In accordance with a preferred embodiment of the invention, a novel

apparatus for and method of simulating physical processes such as fluid flow is

provided. Fluid flow near a boundary or wall of an object is represented by a

collection of vortex sheet layers. The layers are composed of a grid or mesh of

one or more geometrically shaped space filling elements. In the preferred

1 embodiment, the space filling elements take on a triangular shape. An Eulerian

approach is employed for the vortex sheets, where a finite -volume scheme is used

on the prismatic grid formed by the vortex sheet layers. A Lagrangian approach

is employed for the vortical elements (e.g., vortex tubes or filaments) found in

the remainder of the flow domain.

To reduce the computational time, a hairpin removal scheme is employed

to reduce the number of vortex filaments, and a Fast Multipole Method (FMM),

preferably implemented using parallel processing techniques, reduces the

computation of the velocity field.

BRIEF DESCRIPTION OF THE DRAWINGS

The foregoing and other advantages and features of the invention will

become more apparent from the detailed description of the preferred

embodiments of the invention given below with reference to the accompanying

drawings in which:

Fig. 1 is a flow chart illustrating the process flow of a preferred

embodiment of the invention;

Fig. 2 is a flow chart illustrating the process flow of the time marching

step of the preferred embodiment depicted in Fig. 1; Fig. 3 is a block diagram of a computer system in accordance with a

preferred embodiment of the invention;

Fig. 4 illustrates a plurality of vortex sheet layers defined above the surface

of an object in accordance with a preferred embodiment of the invention;

Fig. 5 depicts an adaptively partitioned collection of sub-domains of

source points in accordance with a preferred embodiment of the invention;

Fig. 6 depicts a hierarchical tree representation of the domain depicted in

Fig. 5 in accordance with a preferred embodiment of the invention;

Figs. 7a, 7b, 7c and 7d illustrate multiple levels of an adaptive partitioning

of a domain of source points in accordance with a preferred embodiment of the

invention;

Fig. 8 illustrates the first three levels of a seven level partitioned tree of the

domain depicted in Figs. 7a- 7d in accordance with a preferred embodiment of

the invention; and

Figs. 9a and 9b illustrate a vortex reconnection procedure in accordance

with a preferred embodiment of the invention.

* DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

The invention will be described in detail with reference to the preferred

embodiments illustrated in Figs. l-9b. The invention is described herein in its

preferred application to the simulation of vortex elements in fluid flow about a

solid object. However, the invention may be applicable to any type or

configuration of computing system, method or algorithm used to simulate,

model, analyze, or otherwise predict physical processes that encounter the same

or similar problems overcome by the invention described herein.

A physical process can be accurately and efficiently simulated, modeled,

analyzed, or otherwise predicted in accordance with a preferred embodiment of

the invention by implementing the methodology or process flow depicted in Fig.

1. The process is initiated in step S10 by the input of a set of parameters

indicative of the physical process to be evaluated. (For the purposes of

illustration, the physical process specifically described herein in connection with

the preferred embodiment will be fluid flow.) As will be described below, the

initial parameters may address, for example, the number of vortex sheet layers

defined over the surfaces of the object, the number of field points wit in a sub-

domain in a vortex tube, and other data directed to the resolution, efficiency,

complexity and other criteria of the computation or flow environment necessary

to properly evaluate the subject process. Step S10 preferably includes the process of modeling or defining the

geometry of all the solid surfaces of the object. In accordance with a preferred

embodiment of the invention, the object's geometry is defined as a surface grid

or mesh utilizing one or more different space filling shapes. In the preferred

embodiment, the space filling shape is in the form of a plurality of triangular

elements to represent the solid surfaces in three -dimensions. Preferably, the

surface mesh is laid out, however, one or more refining or improving steps may

be added to increase the accuracy of the surface mesh.

In accordance with a preferred embodiment of the invention, the flow

near the boundary or wall of the object is modeled with a series of vortex sheets.

Preferably, the vortex sheets are organized by two (or any number of) fluid

layers. In the preferred embodiment, a layer grid of the vortex sheets 430 is

constructed by taking the wall-normal lines 420 at the vertices of the surface

triangular elements 400, defining the grid points on these lines 410, and

connecting them so as to generate a prismatic grid, as shown in Fig. 4.

In the preferred embodiment, the layers of vortex sheets are of finite

thickness around the boundary of the object. The precise thickness and the exact

number of layers used are some of the parameters input in step S10, as they tend

to change with the complexity of the object under study or with the desired level

(or order) of accuracy of the computation. For example, when used in the

design of an automobile, the input may specify 10 layers making up the vortex

sheet with a total thickness of 1% of the total length of the automobile.

? Preferably, the layers are of uniform thickness. (In the alternative, the thickness

of the layers may be varied.) In any event, the layers remain fixed overlaying the

surfaces of the object.

Once the geometry and parameter inputs are made in step S10, the

process flow proceeds to step S20 (Fig. 1 ) where initial computations are made .

(and any necessary files restarted) prior to the time simulation of step S30. Thus,

although the vortex method used to simulate turbulence in the fluid flow

marches on in time, there is an initial state of the fluid flow which must be

established in step S20.

The time simulation step S30 represents a plurality of time segments or

steps capturing the change in fluid flow characteristics over time. During each

time step, a sequence of events or sub-steps are performed (S31-S37), as shown

in Fig. 2. In the preferred embodiment, the time steps simulate the creation of

vorticity at the surface of the object, the outgrowth of the vorticity through the

vortex sheet layers, and the release or ejection of the vortices from the outer

sheet layers into the rest of the flow domain. The origin of turbulence and

vorticity, i.e., the instantaneous rotation of the fluid motion, occurs at the

boundary or surface wall of the object. In accordance with the preferred

embodiment, the flow and the creation of new vorticity at the boundary can be

mimicked by stringing together the vortex sheets. A detailed discussion of the

process flow in sub-steps S31 through S37 (Fig. 2) will be described below. After the time simulation step S30, the resulting data produced by the

process flow can be output to one or more applications in step S40. Typical

outputs resulting from the process flow include data on the strength of vortex

sheets, strength and location of vortex filaments, etc. The resulting data can also

be visualized in step S50. Preferably, the visualization is performed using a

three-dimensional visualization tool in accordance with a preferred embodiment

of the invention, as will be described in detail below.

As shown in Fig. 2, the time simulation step S30 is preferably made up of

a plurality of sequential sub-steps S31 through S37. The full sequence of sub-

steps is performed for each time march or step made in the process of simulating

fluid flow using the process flow steps S10 through S50 (Fig. 1), as described

above.

The first sub-step S31 in the sequence for each time step (Fig. 2) of the

simulation involves the computation of the velocity induced on every field point

by every vortex element whose vorticity strength was computed from the

previous time step. The field points are either the nodes of the prismatic layer

grid covering solid surfaces referred to as "grid points," or the nodes of the

vortex filaments referred to as "filament points." The vortex elements are either

the triangular vortex sheets 430 (Fig. 4) that cover solid surfaces 400 or the

vortex filaments. The velocity induced by the triangular vortex sheets on every

field point must first be evaluated. The velocity, u_^ = (u_vn v_vn w_vs), induced by a

to single triangular vortex sheet of constant thickness, d, on a field point is given by

equation (1):

where (x, y, z) are the field point coordinates, (x', y', z') are the

coordinates on the sheet surface, S, and (Ωi, Ω₂, Ω₃) are the vorticity component

distributions of the sheet. The integral of equation (1) is preferably computed in

the following form with equation (2):

after transforming it into the frame of reference where x = y = z = 0 and

the vortex sheet surface is parallel to the xy plane. Note that in this frame of

reference z' < 0 when a field point is above the vortex sheet. The sheets are

assumed to have no vorticity variation across their thickness, while a piece-wise

_linear vorticity distribution is assumed on the sheet surface as follows from

equation (3):

// (3)

Ωⁱ(x', l ) = a_tx' + b_ty' _{+ l}

where a_t, b_n c_I} i = 1, 3 are constants to be determined by the value of

vorticity, Ω at the three triangle vertices, j = 1, 3, where

ttij = a_ix_j' + b_iy_j' + Ci (4)

In this frame of reference, the substitution x' = r cos θ, y'= r sin θ can be

made and the velocity obtained in the form:

/ ..

«>«.

f where

(6) z^z L^l°

(7)

B- r inθdrdθ R + fW+7>\ R (8) h = In sinβdθ o v/r² + z^f' fW+i^a

)

)

For each triangle edge, these integrals are computed after substituting R =

R_mιn/cos (θ _mm - θ) where R_^is the distance of the field point from the triangle

edge. The first integral is computed analytically:

(12)

n The rest of the integrals are computed numerically using integration

techniques such as Simpson, Euler, Trapezoidal rule, etc. In accordance with a

preferred embodiment, a special case is presented. If a field point lies on or very

close to a triangle vertex or edge according to one of the following criteria:

(case 1) R_λ ≤ e;

(case 2) JR.₂ < e;

(case 3) R_z ≤ e;

(case 4) R_mm ≤ 0.05 max {£„ R₂};

(case 5) R_mm < 0.05 max {ϋ₂, R₃}

(case 6) R_mm < 0.05 max {R₃, R,\,

where e- 0.03 R_aVβ and R is the average size of the triangle edges

converging to this node.

In cases 1-3, if the field point is on or near a triangle vertex, preferably

integrals I_x, I Z_^, I and I_xy are not computed on the corresponding triangles

due to the cancellations by the combined effect of all triangles with the common

vertex. If the field point is not on or near a triangle vertex, then z' ≠ 0 and

integrals are computed only along the opposite edge to this vertex.

ιY In cases 4-6, if the field point is on or near a triangle edge, preferably

integrals I_x and I_y are computed without the ln ( \ z' | ) effect due to the

cancellations by the combined effect of the two triangles with the common edge.

If the field point is not on or near a triangle edge, then z' ≠ 0 and, preferably,

integrals are computed only along the opposite two edges to this edge.

The influence of all vortex sheets on a field point is obtained by adding

the effect of all vortex sheets 430 in layers 2 to LYR-1. As will be discussed

below, vortex sheets of layer 1 do not induce any velocity because their influence

is canceled by their mirror-image vortex sheets below the surface (not shown)

that enforce the well-known "no-slip" condition on solid surfaces. Vortex sheets

of layer LTR do not induce any velocity because any vorticity that reaches this

level is transformed into an equivalent filament and its influence is added to the

influence of the other vortex filaments as will be discussed below.

As part of the velocity computation performed in sub-step S31, the

velocity induced by vortex filaments on every field point is computed by first

approximating each vortex tube of each vortex filament by an equivalent point

vortex located in the middle of the tube. The circulation of the point vortex is

equal to the circulation of the vortex tube times the vortex tube length. Then

the velocity on every field point is computed by adding the influence of all vortex

tubes even the ones that were just generated on the upper layer, LYR, of the

grid. In accordance with a preferred embodiment, the summation of all

influences is performed using the well-known Fast Multipole Method (FMM), which is a fast summation scheme for calculating velocities induced by a large

number of vortices. In accordance with a preferred embodiment, the FMM is

implemented using parallel processing to increase the speed and efficiency of

processing, as will be discussed below.

The velocity computations performed utilizing FMM are carried out in

sub-step S32. The FMM is based on combining the velocity fields induced by

collections of nearby vortices into a single local expansion, preferably, in spherical

harmonics. By shifting the center of the expansion to the neighborhood of other

collections of vortices, an efficient velocity field calculation can be effected. To

facilitate the procedure of grouping vortons and field points, the flow domain is

partitioned into cubic boxes. The number and size of these boxes, and,

consequently, the number of vortons per box are parameters that control the

efficiency of the FMM. In a preferred embodiment, the uniform FMM is

employed, where the domain is partitioned into boxes of equal size. As discussed

below, an alternative methodology, the adapative FMM, where the domain is

partitioned into boxes of differing sizes but with a predetermined maximum

number of vortons and field points per box, can also be used. The partitioning

of the flow domain into source boxes S and field boxes F of equal size is achieved

by the following procedure:

• Level 1: The entire flow domain is enclosed inside a cubic box.

• Level 2: The box of level 1 is divided into, preferably, 8 boxes (2 boxes per direction, 2³ = 8) of equal size.

/ι • Level 3: Each box of the previous level is divided into the same number of boxes (e.g., 8 boxes) of equal size.

• Level n: Each box of the previous level is divided into the same number of boxes (e.g., 8 boxes) of equal size.

The FMM reduces the complexity of the O ((N_v)²) algorithm by using

multipole source expansions to compute the influences (field expansions) of the

vortons on the field points, as well as the property that these expansions can be

translated. In the preferred embodiment, multipole source expansions are first

used to represent the influence of vortons within each source box by a truncated

series of spherical harmonics. This multipole source expansion can be thought of

as the equivalent of a single vorton at the center of S whose far field influence is

approximately the same as the influence of each individual vorton in S. A further

reduction in computational complexity is achieved by not applying the source

expansion to each individual field point in F, but instead, on a single point such

as the center of J⁷, for example. This is achieved by shifting the effect of source

boxes onto field boxes via truncated expansions in spherical harmonics. Finally,

the resulting field expansion centered in F can be translated to all the field points

in i⁷, using another truncated series, thus yielding the velocity at each field point.

The truncation order, N, is the same for all these series, and it is the parameter

that controls the accuracy of the FMM in evaluating the velocity field in

comparison to the direct computation.

The above FMM procedure can be used to evaluate the influence of

vortons that are well-separated from the field points. For each field point, the vortons that are well-separated are defined, for example, as the vortons that do

not belong to the same box as the field point or any of the surrounding 26 boxes

that have a common face/edge/vertex.

The parallel implementation of this procedure is based on the following

algorithm:

(a) Distribute boxes and all necessary data (e.g., coordinates of

vortons and field-points within each box) amongst the N_p processors.

(b) On each processor, compute multipole source expansions at boxes 5.

(c) Perform total exchange of multipole source expansion coefficients

between processors.

(d) On each processor -

• Compute field expansions centered in boxes F; and

• Translate field expansions centered in box F to all field points in F. An alternative FMM methodology, referred to as the parallel adaptive fast

multipole method, may be used in accordance with a preferred embodiment. In

the adaptive method, the computational domain is partitioned into a hierarchical

collection of sub-domains which is represented in a tree structure. The shape of

the tree is determined by the distribution of the source points and the constraint

that each of the lowest level sub-domains contain at most a predetermined maximum number of points. When a source point is introduced into the

domain, and the box which contains it already has the maximum allowable

number of points, then the box is sub-divided as, for example, into 2^N equal size

child boxes, where Nis the spatial dimension, and the points in the original box

are re-assigned to the appropriate child boxes. This adaptation gives rise to a

quad-tree in 2 dimensions, as shown in Figs. 5 and 6 and an oct-tree in three

dimensions. With this tree structure, the multipole expansions are used in a

hierarchical way to reduce the complexity of the O ((N_s)²) algorithm to order O

(N_slog ϊNζ), where N_sis the number of source points. This approach is

significantly different than the fast multipole method for the Helmholtz

equation. There, FMM is used to decompose a matrix problem and the resulting

parallel algorithm has different communication patterns and load balancing

issues. This algorithm can be parallelized, for example, by following the Bulk

Synchronous Parallel model for parallel programming as described in "Scalable

Computing," Computer Science Today: Recent Trends in Developments, W.F.

McColl, pages 46-61 (1995), which is herein incorporated by reference.

The adaptive FMM reduces the complexity of the O ((N_s)²) algorithm by

using multipole expansions to compute influences of source points on distant

field points. In a preferred operation, for each leaf box (i.e., the childless boxes

in the tree structure), multipole source expansions are computed. Each

expansion is centered about the center of the box and represents the combined

effects of all source points in that box. The expansions in all the childboxes are

/ f then translated up the tree to the center of their parent boxes. The expansions in

the parent boxes represent the combined effects of all source points contained in

the parent's children. This translation continues up the tree until a certain level

is reached (e.g., level 2; level 0 is the root node). At level 2, all boxes at level 2

and higher (i.e., down the tree) have source expansions which account for all the

source points contained in their portion of the domain. These expansions are

used, when appropriate, to compute the effects of their associated source points

on distant field points.

Preferably, at every field point, there is an accounting of the influence of

all source points. When field points are sufficiently separated from source points,

multipole source expansions are used to represent the combined effects of the

source points. This is preferably done by appropriate translations of the

multipole source expansions to the near-field region of the field points. The

resulting expansions are referred to as the field expansions. The field expansions,

then, represent the influence of distant source points. In general, these field

expansions are translated to the center of field boxes, translated down the tree to

the leaf boxes, and then translated to the field points themselves.

In the preferred embodiment, the field expansion computations start at

level 2. First, field expansions are constructed at the center of all level 2 boxes.

The criteria for which source expansions to use is preferably based on box size

and separation distance, hence level, as will be described below. Next, these field

expansions are translated to the child boxes at level 3. The field expansions at

U> the level 3 boxes are added to by distant source boxes which satisfy the level 3

separation criteria but which have not already been accounted for at a previous

level. This procedure continues down the tree to the leaf boxes. At the leaf

boxes, the field expansions are translated to the field points themselves

accounting for all distant source points. The influence of the near source points

are then calculated.

In the preferred embodiment, this procedure identifies the appropriate

source boxes to use for the multipole field expansions during the downward

traversal of the tree. The union of all source boxes used, including source boxes

of near source points, preferably covers the computational domain without

overlap as each source point preferably is accounted for once and only once. The

tree structure and the rules for applying the multipole expansions are used to

identify these boxes. The rules are preferably based on box size and separation

distance. The optimal choice of boxes will minimize the number of boxes and

maximize the use of the approximations provided by the multipole expansions.

In the preferred embodiment, there are at least three conditions under

which the multipole expansions are used. The conditions are preferably based on

box size and separation distance and are illustrated herein for some distant source

box S relative to field box F.

Under the first condition, 73 is any box, I LB 1 1 can be defined to be equal

to the maximum dimension of any side of B. Where d ( S,F) is the minimum distance between any boundary point of 5 and F, the source box S is well-

separated from field box F if and only if

d(S,F)≥m∞(\\s\\, ILFll).

For example, in Figs.7a-7d, the green source boxes are well-separated

from the gray field boxes. When S is well separated from F, the multipole source

expansion in 5 can be translated to the center of F and eventually used to

approximate the influence of all source points in 5 on all field points in F.

Under the second condition, the separation distance is such that the

multipole source expansions are valid but the translation to the center of the field

box is not. For this case, the source expansion is applied directly to the field

points themselves. Source box 5 and field box .Fare referred to as being source-

separated in this case. Source box S is considered to be source-separated from

the field box F if and only if

IISI IIFII

and

llsll< (5,J < I IF 11.

For example, in the upper left panel of figure 4 the blue source boxes are

source-separated from the gray field box. Under the third condition the separation distance is such that the source

expansions are not valid but the influence of individual source points can be

added to the field expansion at the center of the field box. For this case, source

box S and field box .Fare referred to as being field- separated. Source box 5 is

field-separated from the field box F if and only if

I LFI l lsl l and

\ \F\ \ < d (s, F) < \ \s \ \.

For example, in Fig. 7a the yellow source boxes are field-separated from

the gray field boxes. For each leaf field box, there is an optimal partition based

on these separation rules and the requirement to minimize the number of boxes.

In Fig. 7a, the optimal partition is represented by the union of the gray field box

in the upper left panel, i.e., the generating leaf box, and the red, green, and

yellow source boxes in all panels.

In accordance with a preferred embodiment, a partition of the

computational domain is a set of boxes with non-intersecting interiors whose

closed union is exactly equal to the computational domain. An optimal partition,

from the perspective of the adaptive FMM algorithm, is one which minimizes the

number of boxes and maximizes the use of the multipole expansions. The tree

structure can be used to define many different partitions besides the most natural

one given by the leaf boxes. For example, a leaf box and all its siblings can be

13 replaced by its parent box. If Fis some leaf box containing field points, then

there will be an optimal partition based on separation distance from F which can

be used to evaluate the influence of every source point on the field points in F.

Each leaf box JPin the tree structure generates an optimal partition of the domain

which is needed in the adaptive FMM algorithm to compute influences.

In the described embodiment, the computational domain is decomposed

into non-intersecting open sets such that rules for applying the multipole

expansions apply to each set. The covering of these sets with boxes from the tree

structure identifies the boxes used to apply the multipole expansions. In the

discussion below, the boxes are considered to be open boxes. Closed boxes will

be indicated with an over-bar mark or OB (over-bar) subscript.

The optimal partition of the computational domain generated by a leaf

field box Fis based on the extended colleague lists of and of Fs ancestors. The

extended colleague list of .Fis defined as the list of all boxes adjoining F which

are at the same level as For lower. (The root box is at level 0; hence, lower level

means higher up the tree; hence, a bigger box.) By definition, the boxes in the

extended colleague list completely surround F unless F shares a boundary with

the computational domain. For convenience, F and its ancestors are denoted by

placing the box level in a superscript. For instance, if Fis at level i, then F= F ^l

and Fs parent is F ^ll.

if Where F¹ is the leaf field box generating the partition, {E₁ (F¹)} is the set

of all boxes in the extended colleague list of F and S_OB (F ) is the closed union

of these boxes:

Define a sequence of sub-domains as follows:

τ>^l-^l F^l) = ε(F^{L~ x}) - ε(F^l)

_*>.l^l--2 _/Fτ l^l _\) — ε c ι(F Π^¹--2^'2 N) - ε c /(Fcl^l —^~ 1^l -)

T>°(F^l) = 0.

Since the closure of ε (F ^l) was defined by a collection of closed boxes,

each D (F ^l) is an open set whose interior contains the appropriate boundary

points of the boxes which define it. By construction, the following relations

hold:

5^' D_OB = D (F ^l )

i = 1

D (F ¹) ΓΛ D (F ^x) = 0, i ≠ j

where D is the computational domain.

In the preferred embodiment, this decomposition forms the basis of the

optimal partition of the computational domain with respect to the adaptive

FMM. Preferably, the sub-domain D¹ (F^l) will be partitioned into sets of

source -separated boxes and near-field boxes relative to F¹. The sets D (F^l) for j

< 1 will be partitioned into sets of well-separated and field-separated boxes

relative to F ^,+1.

Where A, (F ^l) is the set of all leaf boxes contiguous with F , ^4_QB (F ^l) is

the closed union of these boxes:

(Fⁱ) = U^^F'⁾-

0 ά Define

C^l = F^d + Λ(F^d).

By construction, L' c D¹ ^ ¹). Also define

The decomposition of D^l (F *) given by

_V ⁱ(F^l) = C^l + C^z

is the desired decomposition for near-field boxes and source-separated

boxes.

The covering of the sub -domain I) with boxes from the tree structure is

unique. Also, because each box in the covering of I) adjoins F ^x, the influence

of all source points in J) on all field points in F ^l must be computed exactly. In

Fig. 7a, the gray box represents F¹ and the red boxes are the unique covering of

L

A covering of the sub-domain If with boxes from the tree, in general, is

not unique. The possibility exists that a box B and all its siblings are contained in

U. In this case, that portion of ³ can be covered by the parent of box B. Using

the lowest level boxes possible to cover If gives the optimal partition. Any box

L 1 in this covering will be source-separated from i⁷ '. In Fig. 7a, the If region is

covered by blue boxes. In this particular example, the covering is unique.

However, it could have been the case where these blue boxes are parent boxes.

The regions I) and If are precisely the regions covered by the list 1 and

list 3 boxes in Strickland and Baty, "A Three-Dimensional Fast Solver for

Arbitrary Vorton Distributions," SAND93-1641 (1993), which is hereby

incorporated by reference. The partition of D¹ (F ^γ) when j < / is similar to the

above. Let

{₍%{ )} = {B \ B is a leaf box, B € {Ei(F and B C V>(F¹)}.

That is, B is a leaf box in the extended colleague list of F which is

contained in the subset D (F ^l). Let C_OB ^J (F ') be the closed union of these

boxes.

Define

L^4J = σ (F

and

£²J= TP(F^l) - C^itj.

r The decomposition of D (F ), j < lis given by

^( ') = £^2,J + ^£4J L^2j n L^4j = 0

is the preferred decomposition for well-separated and field separated

boxes.

The covering of the L^4j sub-domain with boxes from the tree is unique.

Because ^4j is constructed from leaf boxes, they can be used as a covering. If B is

any box in this covering, then it is the case that not all of B's siblings will be in

the covering; hence, no collection of boxes in L^4j can be replaced with their

parent box and so the covering is unique.

If S is any box in the covering of L^4j, then S is a leaf box, and, because it

is also in the extended colleague list of F, it must be the same size or bigger than

box F^J. Because

p³(F^l) [) V^}+l ( ) = ,

L9 S can not be in the extended colleague list of F ^,+1, hence S is separated

from F ^,+1 by a box the size of F ^,+1. Therefore, all boxes covering ^4j are field-

separated from F ^,+1.

In Fig. 7a, L⁴'³ is covered by yellow boxes and is field-separated from the

gray box, F⁴. Similarly, in Fig. 7b the yellow box covers L⁴'² and is field-separated

from the gray box, F³. For the particular choice of F⁴ in Fig. 7c, the L⁴'¹ sub-

domain is empty; hence, no yellow boxes are shown.

The L^2j sub-domain can be covered with boxes the same size as F ^,+1.

Suppose S² L^2j and suppose I |S 1 1 > I IF ^,+1 1 1. Then, if S is a leaf node, S will be in

L⁴'¹ and so cannot be in L^2j. If S is not a leaf node, then descendants of S at level 1

can be used in the covering. Suppose S² L^2j and suppose I |S 1 1 < I |F ^,+1 1 1, then all

siblings of S's must also be contained in L²° and so can be replaced by their

parents.

Any box S in this non-unique covering of L^2j is well-separated from F ^,+1.

This is the optimal covering since both multipole source and field expansions can

be used relative to F ^,+1. In Figs. 7a-7c the green boxes cover L^2j.

The regions L2 and L4 are precisely the regions covered by the list 2 and

list 4 boxes in Strickland and Baty, "A Three-Dimensional Fast Solver for

Arbitrary Vorton Districutions," SAND93-1641 (1993), which is hereby

incorporated by reference.

70 In summary, the partition of the computational domain based on the leaf

field box F 1 is a collection of non-intersecting regions given by:

v = c^l + c³ + (J (c²* + c^AA .

]=2

In accordance with a preferred embodiment, each leaf field box gives rise

to a separate partition. This is the partitioning strategy used in the parallel

adaptive fast multipole code.

The above partitioning strategy was defined from a leaf box up. The

algorithm itself will construct these partitions from the top of the tree down. The

higher up the tree, the more each of the partitions have in common with one

another. For instance, if a sibling of the gray node in Fig. 7a were chosen instead

to generate the partition, the arrangement of the boxes would be different on the

first panel, but not on the subsequent panels since their ancestors are the same.

That is, the partitions given by L²'' and L⁴'^j; j < 4 will be identical; hence, the field

expansions from the top down to level 3 will be identical and only need to be

computed once.

7 / In accordance with a preferred embodiment, the parallel algorithm for the

adaptive fast multipole code may employ a known model such as the Bulk

Synchronous Parallel (BSP) model for parallel programming. The BSP is a

simple programming model based on a superstep structure. A superstep consists

of either calculations or communications, or both; only one sided

communications are allowed. Each superstep is followed by a barrier

synchronization for all processors. The simple structure of the BSP model has

some very important benefits. These include: (1) the task of debugging a

message passing program is much easier; and, (2) one can, with knowledge of the

complexity of the algorithm, predict how a BSP algorithm will perform on a

particular machine without writing a single line of code. The adaptive FMM

code has 2 supersteps: one for computing the source expansions and one for

computing the field expansions and velocities. The parallel algorithm, including

the details of these 2 supersteps, are described below.

When solving most n-body problems, the computational costs far exceed

the memory requirements. This is true for the vortex method, even for O(10⁷)

source points. Because of this, and because of the complexity of parallelizing the

adaptive FMM code, certain simplifying assumptions are made in this illustrated

embodiment. These assumptions are: (1) each processor has its own copy of the

full tree; (2) a total exchange of the source expansion coefficients is done at the

conclusion of the source expansion super-step. (Preferably, the field expansion

coefficients may remain local to each processor to remove the need to exchange this data.) These assumptions make it easy to distribute work to the processors -

only node numbers of starting nodes need to be passed to each processor to

indicate its responsibility.

The parallel algorithm for computing the source expansions is:

(a) Partition the tree.

(b) Assign subtrees to processors.

(c) Source Expansion Superstep: On each processor

(1) For each sub-tree

• compute multipole source expansions at all leaf boxes; and

• translate the source expansions up to the root of the sub-tree (level l_r(T,)) or to level 2 as appropriate.

(2) Do total exchange of multipole source expansion coefficients.

(d) If necessary, translate source expansions from level l_r(T,) to level 2.

Step (a) partitions the tree into a collection of sub-trees for distribution

across the processors. The main objective of this step is to load balance. The

original tree is pruned in such a way that the resulting sub-trees have as close as

possible the same number of source points. In figure 5, the first 3 levels of a 7

level tree are displayed. The tree is partition into 22 sub-trees. The root node of

13 each sub-tree is colored blue. The sub-trees are distributed to the processors in

step (b) and the source expansion superstep is done in in step (c). The special

sub -tree consisting of the original root node and extending down to but not

including the root nodes of all other sub-trees is referred to as the root-node

sub-tree. In Step (d), all processors operate on this sub-tree which results in a

small amount of redundant calculations.

The number of sub-trees, T„ is a function of the number of processors,

the distribution of the source points, and the desired granularity. For higher

granularity, many more sub-trees than processors are needed. The advantage of

high granularity is load balancing is easier to achieve; however, more sub-trees

means the root node of each sub-tree is further down the original tree which

requires more redundant calculations at the root-node sub-tree. Once the sub¬

trees are identified, they can be distributed to the processors in a number of

different ways: (1) a contiguous allocation based on the number of sub-trees and

the number of processors, (2) a round-robin distribution, and (3) a more

sophisticated technique based on solving the knapsack problem. These

characteristics of the parallel algorithm can be controlled with inputs to the code.

The parallel adaptive algorithm for computing the field expansions is:

(e) Field Expansion Superstep: On each processor

-i (1) For level 1 = 2 to the level just above the leaf nodes, do for each sub¬

tree

• for each non-leaf field point box F at the current level; do

- construct list 2 and list 4 boxes relative to field box F¹

- compute field expansions centered in box F¹ from the list 2 and list 4

boxes

- translate field expansions to the children of box F¹.

(2) For each leaf field point box F¹; do

• translate field expansions centered in box F¹ to all field points in Fl

• construct lists 1 and 3 boxes for box F¹

• add influences of all source points in the list 1 boxes and box F¹ itself to each field point in F • add influences of the source expansions in the list

3 boxes to each field point in F¹

• compute velocities at the field points.

(3) Do a total exchange of velocities.

In the preferred embodiment, the algorithm for computing the field

expansions starts level 2. Depending on a number of factors, there may be

? r insufficient parallel work at this level and so some redundant calculations may

have to be done. This depends on the largest level of the root node of the all the

sub -trees; refer to this level as l_max. All processors start the computation at level 2

and proceed down the tree computing all field expansions until level l_max is

reached. From this point forward, only the nodes in the sub-tree are processed. If

L_a > 3, then redundant calculations at the top of the tree are done. If l_max < 2,

then no redundant calculations are done.

The computation of the field expansions can be accomplished in one

superstep. The field expansion data at the leaf boxes do not need to be

communicated to their neighbors. It will be used only locally. The superstep

ends with a total exchange of the velocities.

In the preferred embodiment, the influences of all the source points in D

are accounted for only once. The field expansion calculations begin at level 2 and

continue down the tree. In Fig. 7a, this can be seen graphically for a particular

field box F¹ by the fact that the entire domain is ultimately covered with red,

green, and blue boxes, plus F¹ itself. Step (1) of (e) computes the field expansions

based on the green and yellow boxes while step (2) computes the influences

based on the red and blue boxes, and based on the source points residing in the

leaf field box itself.

3 N_υ(S) Number of source points in box S s Cartesian coordinates for the center of box S x; = (xi, yι, z Cartesian coordinates for the ^tth source point in S

G(S); Vector valued stength of the ith source point in 5^" p( ; — s-), α_i(x_i — s), β (x.ι — s) Spherical coordinates of ith source point in S with origin at the center of S y^m(α(x_j - s), 3(x; - s)) Spherical harmonic

Vector valued coefficients of the expansion centered at s

Field point

P

Order of the expansion

N

Source Expansion Coefficients

Expansion coefficients at the center of source box S due to individual

source points in S are:

N,(S)

A(s)J = ∑ G,(S xι - s)J Y_j-^,'(α(xι - s),/?<Xi - s))

1=0

Lemma: _Ak _ _A-k _{where A}-k _{is the} complex conjugate of A_j ^k.

The lemma follows immediately from the fact that G, and p_jA are r< 17 =^y,^"

Translation of source domain expansion coefficients centered at box S to the

center of box P are computed by

3 1 where ml = _maχ(-„, _K _ _N^ _{and mh =} min(n, JK -f- jy) _t

Lemma: A* = A^

From the definitions of D_m/n and Ja_m/n , we have the following:

—m

K = On

— m

J % = JaZ^™

Using this and the fact that p ₌ A -: k , it is apparent that for some

integer K such that 0 < K < N ,

1 T

J m=-K+N τ_n-K-mr₎mr₎-K ri-K- y-im+K)/ /\ n

Σ ' ^{m u}n^uj J-n 3+n ' / P _ΛK+Γ

2^ f)k ^Aj-n n=0 τn=— n J

= ^J" ^m Ϋ⁼"ⁿ Ta ^{mK~m u}D^mυn3^KU n_d ^K-^~n^mY ^l j^m+^~n^K (a^a^P β)⁾ P oⁿ βK-txx _ K n=0m=K-N 3

A similar argument applies for negative K.

Because of the last 2 lemmas, there is only need to save the expansion coefficients with positive k indices.

Field Expansion Coefficients

The local expansion centered at some field point p due to source points in a well separated source box S is given by

?f

(_l)n(_l)mm(H,|fc|)₎ if m • fc > 0 (-1)" otherwise

and due to N_v single source points from source box S is

The translation of the field point expansions from point p to point q is:

where

(_l)ⁿ(-l)^m, if » /c < 0

Jb. n,m (-l)ⁿ(-l)^fe_m if m * k > 0 and |fc| < |τn| (— l)ⁿ otherwise

Lemma: B^ = B,^~k

Lemma: B^ = B^~A

In accordance with a preferred embodiment of the invention, the

strength of the triangular surface source panels can be computed by the

numerical solution of the following integral equation, which represents the no-

penetration condition on solid surfaces,

T where the left hand side is the surface normal velocity induced on the

surface point (x, y, z) by a surface source strength distribution q(x', y¹, z'), and

the right hand side is the opposite of the surface normal velocity induced by all

the vortex elements of the flow. In this equation, r = (x - x¹, y - y', z-z' ) and n is

the surface normal unit vector.

To facilitate the numerical solution of equation (13), the no-penetration

condition is enforced at the node points of the triangularized surface and results

in the following linear system of equations for the source strengths at the same

node points:

A, q, « - (n • u,), (14)

where each element of the array A,, is equal to the normal velocity induced

at surface node point i by a piece -wise linear source distribution with unit

strength at surface node point j and zero on all other surface node points. The

construction of matrix A^ is performed by looping over all surface triangles

assuming a piece -wise linear source distribution where, in turn, each triangle

vertex has unit source strength while the other two have zero, and computing

the normal induced velocity on all the surface node points.

The construction and inversion of array A, may be performed in the

precomputation step S20 (Fig. 1) of the process. During each time step, a new

i 01/01237

panel strengths that enforce the surface no-pene.ation conditⁱon.

Once the source strengths are computed in sub-step S33, the velocⁱty, u_s

veloc_ity u. induced by all vortex elements (as computed after stages 1 and 2) to

result in the total velocity field in sub-step S34.

,_{n v} w ) _induced by a single triangular surface source The velocity, u_s = (u_s, v„ w , mu^{u j}

panel on a field point is given by the equation.

where (x, y, ) are the field porn, coordinates, <*' , /, -') « the _{l cr}fcre S and Q is the source strength coordinates on the source panel surface, S, ana ^ ral of equation (15) is computed in the distribution of the panel. The integ

following form:

the frame of reference where x = y = z = 0 and after transforming it into _Λe source panel surface ,s p r* ro the xy plane, and assum^tng a

the panel surface. The constants a, b, are linear source strength distribution on

3 determined by the value of source strength, Q_j7 at the three triangle vertices, j =

1, 3, where

(17)

Q_j = ax_j' + by_j' + c

As described above, in this frame of reference, the substitution x, = r cos

θ, y, = r sin θ can be made, and the velocity obtained in the form:

where the integrals are as defined previously above. For efficiency reasons,

the computation of induced velocities as described above is carried out for each

field point and the computed velocities are then interpolated to calculate the

velocity induced on filament end points which lie inside the vorticity layers.

Filament end points which lie outside the vorticity layer are assigned approximate

velocity field which is computed assuming that the vorticity of each surface

triangle is concentrated at its center. Each time step, preferably an oct-tree-based

search is executed to identify the filament points which lie inside the layer and

the triangles which need to be used for interpolation.

Yf A new value for the surface vorticity is computed by enforcing the no-slip

condition on the surface node points. This is achieved by considering an explicit

finite-difference approach in the following form:

ω -f l __ n x u (19)

Ah

where w is the surface vorticity (first grid layer), n is the unit surface

normal vector, u is the velocity induced at the second grid layer and Δ h is the

second layer thickness.

Once vorticity is generated on solid surfaces, its evolution on the layer

grid points is followed by solution of the vorticity equation

dω l (20)

Ϊ - ^{V x f +} B^Λ

where F = u x w is helicity. In the prefered embodiment, an explicit

scheme may be used for the solution of the vorticity equation. Thus, in sub-step

S35, the derivatives of helicity on every grid point are computed using the

divergence theorem on the surrounding prismatic volumes

(21)

— dF_r 1 J f_sF_injdS where S is the total surface and V is the total volume. A similar approach

can be used for the diffusion term, although for high Reynolds number flows

only the wall normal diffusion term is retained and computed by a finite-

difference scheme.

In accordance with a preferred embodiment of the invention, new vortex

tubes are created in sub-step S36 out of the vorticity residing in the outermost

layer of the triangular prisms in the mesh surrounding solid bodies (Fig. 4). The

algorithm precedes as follows: Assume that the vorticity in the outer layer of the

prisms has been set to zero at the beginning of a time step. Then, after solution

of the finite volume approximation to the vorticity equation, new vorticity

appears in the top layer of triangular prisms, say Ω, in the ith prism. If the

vorticity held in the ith prism exceeds a critical value, say "vrt-thrhld," i.e.,

|Ω, I > vrt-thrhld, then either this vorticity is used to create a new vorton at this

location, or else it is added in to a pre-existing vorton which was created at this

location during an earlier time step. If the vorticity magnitude falls below the

threshold no contribution to the vorton population is made from this location

during the time step.

A completely new vorton may be added if either there has not recently

been a pre-existing vorton at this location, i.e., during the previous time step |Ω,

was below the threshold, or else a pre-existing vorton at this location has just

been "released," i.e., set free to join the general population of vortons so that it

f is no longer associated with its birth location. Pre-existing vortons are "released"

for at least three different reasons:

1. their strength exceeds a criterion "max-crcltn,"

2. they have convected too far away from the center of the prism from

which they originated (the distance is set to be equal to the average size of the

sides of the triangle associated with the prism), or

3. |Ω, I drops below the threshold for new vortons during a time step.

If a new vorton is to be added, it is preferable that it be specified at least

according to its position in space, its strength, i.e., circulation, its orientation and

its length. If the vorticity at a location is to be added to a pre-existing vorton,

then the position, strength, orientation and length should be specified for the

combined vorton.

New vortons are introduced at the center of the ith triangular prism. They

are oriented in the direction of the vorticity field they are to represent. If Δ s

represents the length of the new vorton, T denotes its circulation and V_;

represents the volume of the ith triangular prism, then we must have

(22)

Ωi [ Vi = FAs.

This guarantees that the triangular prism and the vorton replacing it, will

have the same far field velocities.

Y? According to equation (22), only the product of vorton length and

circulation must be specified, not eachindividually. In accordance with a

preferred embodiment, a procedure has been developed to first select Δ s and

then pick T so as to satisfy equation (22). The vorton length is determined as the

minimum of either the average length of the sides of the triangle associated with

the z^'th prism, or the distance along the line between the top and bottom surface

of the triangular prism where they are intersected by a line through the center of

the prism in the direction of the vorticity.

If the new vorticity is to be added to a pre-existing vorton, (i.e., a vorton

that has been created earlier at this location but has not yet satisfied the release

criteria), then the following procedure may be followed: The new circulation of

the vorton is taken to be the sum of the pre-existing circulation plus the new

circulation determined from (). The new position of the endpoints of the vorton

are found by taking a weighted average of the end positions of the new and pre-

existing vortons, with the weighting done according to the circulation of the new

and pre-existing vortons. Specifically, (x,; y_x; z_x) and (x₂; y₂; z₂) represent the

beginning and end points, respectively, of the pre-existing vorton, and (a,; b,; c. )

and (a₂; b₂; c₂) the similar points for the new vorton. Then the position of the

pre-existing vorton is modified as follows:

(1 - a)yj + abj

yr where j = 1, 2 and α = r„e_W/rβW.

At every time step, the velocities at the locations of the end points of all

vortons are calculated, and new updated positions of the vortons are found via a

simple Euler scheme in sub-step S37. Specifically, if xn represents the endpoint of

one of the vortons at the start of the nth time step of the calculation, then

_χn+l _ _χn _{+ u}n_At (24)

gives the new position of the vorton end point. Here u" is the velocity

field calculated at the point xⁿ at the nth time step. Alternative schemes for the

time advancement of the vorton end points may be implemented, for example,

schemes with greater accuracy. As the result of the convection motions

described above, vortons will tend to change in length. If an individual vorton

grows in length beyond a certain threshold, then it is divided into two separate

vortons of equal length. The two new vortons share a common point at the

computed midpoint of the two end points of the original vorton. If a particular

vorton reduces in length to a size falling below a threshold then it is combined

with one or the other adjacent vortons in the filament.

ff To model viscous diffusion as it effects vortex tubes in the calculation, the

following method has been implemented in accordance with a preferred

embodiment: Let r_c denote the equilibrium radius for vortex tubes with a

Gaussian core in which the effects of vortex stretching are in equilibrium with

viscous diffusion. For a tube whose length increases from Δ s to Δ s' in time Δ t

it may be shown that

where R,. is the Reynolds number.

For each individual vortex tube at each time step, an estimate of the local

stretching ratio Δ s'/Δ s is made by averaging over the stretching felt by each of

the vortons composing the vortex tube. This average is used to get a value of r_e

pertaining to the vortex filament. The actual radii of the vortex filaments used in

Vorcat, r, are the same constant value for all vortices. For any one tube the

computed value of r_e may be larger or smaller than r.. If r. > r_e, for the given

amount of stretching experienced by the vortex tube, then its equilibrium radius

is smaller than the actual radius it has in the computation. Hence, if given the

opportunity, the vorton will tend to get thinner. On the other hand, if r, < r_e, the

vortex is not experiencing enough stretching to justify its radius, and if given the

opportunity it will diffuse outward and become larger.

f? In accordance with the preferred embodiment of the invention, these

differing conditions are modeled by the following mechanism: If r, < r_c, a loss of

circulation from a tube is implied. The amount lost is determined by the extent

to which vorticity would diffuse beyond r. due to the imbalance of vortex

stretching and diffusion. The necessary relation is

_rn_{+ 1 = r l} 4Δt ^ 1 - - 1 _{) }} (26)

where Tⁿ is the circulation at time n. If r. > r_e, the circulation is left the

same, i.e.,

pn-h l _ γ,n _{^ (27)}

and if the tube faces a net contraction, i.e. if Δ s' < Δ s, then the

maximum circulation loss is assessed, namely

-.n+l ^Δt2Λ (28)

As vortex filaments and their component vortons evolve during the

calculation, at least two techniques to limit the development of small scale flow

features may be used. These are the processes of hairpin removal and vortex

reconnection.

Where s, denotes the vector connecting the end points of the ith vorton,

and two adjacent vortons s, and s₂ are lying on the same filament, the angle θ

formed between these adjacent vortons can be obtained from the identity Si ^• S₂

COS 0 =

Isil I (29)

When θ is in the vicinity of π, the filament has a discernible kink, fold or

"hairpin" at the point between the two vortons. In accordance with a preferred

embodiment, such "hairpin" configurations are removed from the filaments

when cos θ falls in magnitude below a threshold, e.g., cos θ^* where θ^* is an angle

near π . The hairpin is removed by replacing the two vortons by a single straight

vorton connecting the outer end-points of the original vortons.

Vortex "reconnection" pertains to the situation when two nearby vortons,

for example, say s, and s₂, belonging to different filaments with very similar

circulations, have nearly equal and opposite rotation, i.e., s, is approximately anti-

parallel to s₂. In this case, the rearrangement of vortons and filaments as shown in

Figures 9a and 9b. This procedure will allow a single vortex ring to split into

two separate rings, or allow two separate rings become one ring. In accordance

with a preferred embodiment of the invention, vortex filaments are not

constrained to be in the form of rings, so vortex reconnection acts to cause the

merger or separation of individual filaments.

Similarly, in accordance with a preferred embodiment, vortex

reconnection may also be implemented as follows: The circulations of two vortons from different filaments are deemed to be within a certain percentage of

each other (e.g., 10%). Then, letting x, and x₂ denote the center points of the

vortons s, and s₂, respectively, x₁₂ = x₂ - x, is the vector connecting the centers of

the two vortons. An estimate of the distance, (e.g., d) between the vortons is

made which allows the decision to be made as to whether vortex reconnection

between these two vortons should be pursued.

d is determined as the average distance:

where d, is the length of the vector x₁₂ projected onto the direction of a

vector perpendicular to s₂ and passing through x, , while d₂ is the length of the

vector x₁₂ projected onto the direction of a vector perpendicular to s, and passing

through x₂.

When d is less than a threshold value, e.g., half the radius of a vorton,

then the angle between the two vortons is computed using equation (29), as in

the case of hairpin removal. If this angle is close enough to ss then the vortons

are considered to be anti-parallel, and then the vortex reconnection shown in the

figures is implemented.

Γ In accordance with a preferred embodiment, at least two different

methods can be implemented for modeling the process by which vortices reach

the end of their useful existence in the calculation:

1. During each time step of the calculation, if, as a result of subdivisions

taking place since its creation, the number of vortons composing a particular

tube exceeds a critical number, then that vortex tube and all of its component

vortons are eliminated from the computation.

2. During each time step of the calculation, if, as a result of the diffusion

model described above, the circulation of a tube falls below a specified value,

then that vortex tube and all of its component vortons are eliminated from the

computation.

Vortex tubes that have been released from vortex sheets sometimes move

back into the vortex sheet layer. For these tubes, the vortex method described

above may not accurately evaluate the velocity. As a result, in accordance with a

preferred embodiment, the velocity of the re-entered tube is determined through

interpolation from the surrounding triangular vertices. As the velocities on the

edges of the prismatic grid in which the re-entered tube is located are already

Icnown, this interpolation can yield a very accurate reading of the velocity of the

re-entered vortex tube.

ST In order to accurately evaluate the velocity, however, the re-entered

vortex tube must be precisely located within the vortex sheet layers.

Accordingly, a search mechanism such as the well-known oct-tree is utilized. In

the preferred embodiment, the oct-tree is used to search the layered sheets to

locate the re-entered tube. Once the prism that contains the re-entered tube is

found, the interpolation calculation can be performed.

The process flow depicted in Fig. 1 may be implemented in computer

hardware, software, or both. As shown in Fig. 3, an exemplary computer system

is provided having a plurality of modules for performing the various functions

depicted and described above in connection with Figs. 1 and 2. The computer

system, for example, includes input device 500 for entering the geometry and

parameter input information (e.g., Reynolds number, angle of attack, preferred

number of layers in vortex sheet, etc.) needed to ensure the accurate and efficient

operation of the simulation. The input device may be an input file reader,

graphical user interface, keyboard, or any other known mechanism for inputting

data or information. Among other things, the input information is used by

surface grid generator 510 and vortex sheet construction unit 520 to model the

geometric shape of the object and develop its prismatic grid.

Vortex tube generator 530 models the vorticity created at the surface of

the object and tracks its growth outward through the different vortex sheet layers

θ^' until it becomes a vortex tube released from the upper vortex sheet layer into the

remaining fluid flow. The vortex sheet velocity processor 540, FMM processor

550 and source panel velocity processor 560 operate to compute the velocity and

vorticity of field points throughout the flow domain taking into account the

various influences described above. A filament destruction/removal processor

570 is provided to eliminate wasteful computations (e.g., performing hairpin

removal, vortex destruction, etc.) made on vortices not having significant impact

on the calculations of the field point velocities, as described above. Output and

visualization tools 580, 590, respectively, are provided to facilitate use of the

resulting data produced by the simulation of fluid flow in accordance with a

preferred embodiment of the invention.

Visualization tool 590 is preferably embodied as software forming a part

of (or operating in conjunction with) the turbulence simulation computer

system. Visualization tool 590 may also be implemented, however, as part of a

separate computer system or computer hardware. To the extent needed, the tool

will receive from the simulation computer system (or independently calculate)

data structures needed to display the simulation results. For example, the tool

utilizes certain configuration data regarding the underlying geometry of the

object and the vortex sheets (e.g., mappings of the object geometry to associated

triangles; computed normal, velocity, and vorticity of each vertex; calculation of

surface normals; etc.). Typically, the collection of vortons simulated are

encapsulated into a layer. The visualization tool 590 utilizes data on the velocity and vorticity at the two vertices of each vorton in the vorton layer, and

determines if penetration into any of the solid surfaces of the object has been

made by the vortons.

Visualization tool 590 preferably renders data from the simulation in

three-dimensions using different color (or gray-scaling) schemes. The surface

primitives (e.g., space-filling objects such as triangles) have vertices that are

distinguishable in terms of velocity, vorticity, etc. by color or gray-scale.

Aesthetic parameters such as shininess and ambience may be preset or user

configurable. In the preferred embodiment, the vortons are rendered as lighted

geometrical constructs (e.g., cylinders). The size and shape may be scaled in

proportion to given data output from the simulation. For example, the radius of

a rendered cylinder may be scaled according to a representative body size. The

constructs are preferably organized in an octree (or other) data structure to

improve rendering performance. Various parameters of the constructs such as

spin and magnitude may be visualized with use of color (or gray-scale) schemes.

As another example, parameters such as velocity and vorticity may be represented

as scaled, normalized vectors with cone-shaped arrowheads.

In a preferred embodiment, the user interface for the visualization tool

590 is Tcl/Tk-based with pull-down menus to toggle visualization options.

Each layer, along with velocity and vorticity, for example, can be toggled on/off.

The interface permits selective rendering through input of threshold values used

to filter out data not meeting the input thresholds, or through the selection of elements on the display. For example, the user may wish to visualize only those

vortons whose strength exceeds a given value, or the user may simply select or

highlight vorton or body layers to receive detailed information regarding the

selected layers.

It should be readily understood that the components depicted in Fig. 3

may be fully or partially combined into one or more different processing units.

The functionality of each component may be easily implemented separately or in

various component combinations in one or more general purpose computers

executing specific computer software or algorithms. Moreover, the invention

may be embodied as an article of manufacture in one or more computer

programs stored on one or more recording mediums (e.g., floppy or optical disk,

read access memory (RAM), read only memory (ROM), disk drive, etc.).

While the invention has been described in detail in connection with the

best mode of the invention currently known, it should be readily understood that

the invention is not limited to the specified embodiments described herein.

Rather, the invention can be modified to incorporate any number of variations,

alterations, substitutions or equivalent arrangements not heretofore described,

which are commensurate with the spirit and scope of the invention. For

example, although no specific object was described herein as serving as the

surface structure upon which turbulence was simulated, it should be recognized

that any size, shape, or construct would suffice such as a prolate spheroid,

cylinder, etc. The object used in the preferred embodiment is merely a generic

f* representation of the infinite variety of structures that may be the subject of

analysis (e.g., airplane wings, rotors, cabins, automobile exteriors, interiors,

engine components, etc.). The fluid disclosed may be liquid, gas, chemical

mixture, energy, or any process flow subject to evaluation.

The computer system described herein may be centrally located, or

remotely distributed over a network of processing units. The computer system

described in the preferred embodiment is merely representative of the multitude

of systems which may find use for the invention. For example, the computer

system may be part of an automobile design system used to analyze external

aerodynamics of drag forces, noise generation, cross wind stability, etc., during

development of new automotive shapes.

Although the detailed description herein has focused on the simulation of

turbulence using three-dimensional modeling, it should be readily apparent to

those skilled in the art that the invention may similarly be embodied in a system

or algorithm utilizing only two-dimensions.

f/

Claims

What is claimed is:

1. An apparatus for visually simulating on a display device, during

successive time steps, a turbulent flow of a fluid relative to an object, the

apparatus comprising:

modeling means for defining a model of the surfaces of the object;

layering means for defining a plurality of vortex sheet layers surrounding

the surfaces of the object

vortex tube tracking means for defining a plurality of vortex tubes

growing out of an outermost layer of said plurality of vortex sheet layers relative

to the surfaces of the object; and

visualization means for visualizing the object, the layered vortex sheets,

and the vortex tubes during successive time steps.

2. The apparatus for visually simulating on a display device as recited

in claim 1, wherein said modeling means defines a three-dimensional model of

solid surfaces of the object; and

wherein said layering means further comprises a surface grid generator

means for forming the vortex sheet layers by defining a mesh of triangular-

shaped space filling elements stacked over each other beginning with a surface

triangle layer overlaying the surfaces of the object.

ά «

3. The apparatus for visually simulating on a display device as recited

in claim 1, further comprising:

grid means for constructing a layer grid of the plurality of vortex sheet

layers by taking wall-normal lines at the vertices of the surface triangles, defining

grid points on the wall-normal lines, and connecting the grid points so as to

generate a prismatic grid;

vortex tube tracking means for defining a plurality of vortex tubes

growing out of an outermost layer of said plurality of vortex sheet layers relative

to the surfaces of the object;

first field point calculating means for calculating velocity induced by each

vortex sheet on all of the grid points and all points making up said plurality of

vortex tubes;

second field point calculating means for calculating velocity induced by

each vortex tube on all of the grid points and all of points on said plurality of

vortex tubes; and

summation means for summing all resulting calculations from said first

and second field point calculating means using a Fast Multipole Method (FMM)

algorithm, where said FMM algorithm is implemented using parallel processing.

4. The apparatus of claim 3, wherein said summation means takes the

resulting calculations made by said first field point calculating means for all the

vortex sheet layers overlaying the surface triangle layer.

/

5. The apparatus of claim 3, wherein said layering means further

defines a subsurface vortex sheet layer below the surfaces of the object so as to

enforce a no-slip condition on the solid surfaces of the object.

6. The apparatus of claim 3, wherein said summation means employs

an adaptive FMM algorithm in which each vortex to be evaluated is partitioned

into a hierarchical collection of sub-domains represented in a tree structure, for

each field point introduced into the tree.

7. The apparatus of claim 1, wherein said visualization means

provides three-dimensional rendering of selected vortical elements, wherein

detailed information is provided regarding only the vortical elements selected.

67-