Apparatus and Method for Predicting Flow Characteristics

(1) Field of the Invention

This invention generally relates to the analysis of fluid flow past an object and more particularly to a method and apparatus for predicting the characteristics of a fluid flowing past such an object.

(2) Description of the Prior Art

Understanding the characteristics of fluid as it flows past an object, such as an airfoil, is important both from the standpoint of understanding and improving the designs of such objects and in understanding the nature of any turbulence introduced as a result of relative motion of a fluid an airfoil, either by moving of the airfoil through the fluid or by moving the fluid past the airfoil.

In the past understandings of fluid flow have been derived from the observation of fluid flow past a model and by specific measurements. For example, U.S. Pat. No. 3,787,874 to Urban discloses a method for making boundary layer flow conditions visible by applying to the surface of a moving or stationary structural body to be exposed to the flow a reactive layer of at least one chemical color indicator, such as a cholesterinic liquid. The body is exposed to a flow of gas, such as air, which contains a reagent. The chemical color indicator can also be applied together with gelling means and a moisture binder. The chemical color indicator can also be absorbed by a high-contrast, absorbent paper which is then applied to the body. A metal or plastic foil coated with a binder and/or indicator can also be used for this purpose. A boundary layer flow pattern image is produced, which can subsequently be recorded.

U.S. Pat. No. 3,890,835 to Dotzer et al. discloses another approach to chemically recording flow patterns by treating the surface to form a reactive layer, entraining in the fluid a reagent compound which is capable of chemically changing the reactive layer, and then passing the fluid over the reactive layer which is to be examined. In this particular disclosure, a blade or other member of aluminum to be examined is treated to form a thin oxide film by anodic treatment. This film is impregnated with an organic dye. As an air stream containing a reactive substance passes over the treated blade, the acid vapors react with the dye and/or the oxide layer and produce a visible pattern upon the blade. This pattern is characteristic of the boundary layer flow of the air stream. An examination of the visible pattern helps to determine the proper design and operating characteristics of the blade.

U.S. Pat. No. 4,380,170 to Dotzer et al. discloses another process for chemically plotting the boundary layer flows over uncompacted, coated, anodically oxidized aluminum surfaces by using a colored or uncolored liquid or a coating or pointillization with a substance, preferably a dye, soluble in water or organic media and which can be included or adsorbed in the eloxal layer.

U.S. Pat. No. 4,727,751 discloses a mechanical sensor for determining cross flow vorticity characteristics. This sensor comprises cross flow sensors which are non-invasively adhered to a swept wing laminar surface either singularly, in multi-element strips, in polar patterns or orthogonal patterns. These cross flow sensors comprise hot-film sensor elements which operate as a constant temperature anemometer circuit to detect heat transfer rate changes. Accordingly, crossflow vorticity characteristics are determined via cross-correlation. In addition, the crossflow sensors have a thickness which does not exceed a minimum value in order to avoid contamination of downstream crossflow sensors.

These prior art approaches present visualizations or measurements that define certain aspects of the characteristics of fluid flow. However, they are designed primarily to determine characteristics at a boundary layer or some other localized site. Each requires the production of a physical model and physical testing of such models. Moreover, if the testing suggests any change to the shape of an airfoil, it is generally necessary to modify the physical model and run the tests again in order to validate any change. Such testing can become time-consuming and expensive to perform.

More recently, it has been proposed to utilize computer modeling techniques to produce such fluid flow analyses. Such computer modeling is attractive because it eliminates the need for physical models and holds the opportunity to reduce testing, particularly if design changes are made to an object undergoing test. Initially such techniques were applied to circular cylinders using a small number of discrete point vortices.

Eventually additional studies determined that vorticity was useful as a basis for understanding fluid flow. Vorticity is produced at a solid boundary because at the surface the fluid has no velocity (i.e., the fluid exhibits a no-slip condition). Once generated at the surface, vorticity diffuses into the volume of the fluid where it is advected by local flow. Conventional vortex methods generally mime this process. In accordance with such methods, the strengths of the vortex elements or segments originating on the body surface are determined by requiring that the velocity induced by all the vortex elements on the surface be equal and opposite to the velocity at the surface. It is assumed that this vorticity is contained in an infinitely thin sheet at the surface. In these methods a resulting matrix equation is solved for the surface vorticity at all points on the body simultaneously. Vorticity transfer to the flow is then accomplished by placing the vortex elements above the surface.

It has been recognized that these vortex methods have several shortcomings. When computational methods use point vortices in their simulations, mathematical singularities can produce divergent solutions. This has been overcome by using a kernel function that contains a regularized singularity. However, this kernel function depends in certain ad hoc assumptions such as the value of the cutoff velocity and core radius. While the no-slip and no-flow boundary conditions provide information regarding the strength of the surface vorticity and subsequent strength of the vortex element, their use often neglects the effects of all other vortex sheets on the surface. Other implementations of such methods neglect the effects of coupling between the surface vortex sheets and surface sources. Finally, many methods assume a priori a separation point to analyze shedding of vorticity from the surface into the flow that generally requires experimental knowledge of the flow.

More recent prior art has utilized computer modeling based upon the nature of vortex elements at the surface of an object, such as an airfoil. These models then track the motion of each element as it moves into the flow over time to calculate the velocity of each element. While this prior art produces acceptable results, the direct calculation of the velocity of each vortex element produces an N.sup.2 increase in the required time for processing where N is the number of vortex elements for each time step. Such increases can become unacceptable when high resolution demands the calculation of a large number of vortex elements.

Therefore, it is an object of this invention to provide an improved method and apparatus for predicting the flow of fluid past an object.

Another object of this invention is to provide an improved method and apparatus for predicting the flow of fluid past an object that is readily adapted for computer simulation.

Still another object of this invention is to provide an improved method and apparatus for predicting the flow of fluid past an object that minimizes the assumptions used in the predictions.

Yet another object of this invention is to provide an improved method and apparatus for predicting the flow of fluid past an object for a large number of points in an area of interest thereby to provide maximum resolution for the prediction.

In accordance with a method and apparatus of this invention fluid flow characteristics are predicted by defining a model of the object and defining a plurality of vorticity segments over the area of interest including vorticity segments at the surface of the object having a known velocity. In addition, there are defined over the area of interest a plurality of sets of boxes. Each set has a predetermined relationship in size and position to the boxes in the other sets. The velocity for each vorticity segment is calculated based upon the vorticity of that segment, the vorticity of each segment in the same box and in neighboring boxes and the vorticity of groups of other boxes surrounding the neighboring boxes taken collectively. Each of these velocities is summed to determine the velocity of each vorticity segment based upon the influence of all the vorticity segments in the area of interest and to provide information for displaying a representation of the fluid flow characteristics over the area of interest.

The appended claims particularly point out and distinctly claim the subject matter of this invention. The various objects, advantages and novel features of this invention will be more fully apparent from a reading of the following detailed description in conjunction with the accompanying drawings in which like reference numerals refer to like parts, and in which:

FIG. 1 represents a line segment representation of an airfoil cross-section;

FIG. 2 depicts the geometry of a vorticity segment useful in understanding this invention;

FIGS. 3A through 3D depict a box numbering scheme that is useful in implementing this invention;

FIGS. 4 depicts, in schematic form, an apparatus embodiment of this invention; and

FIGS. 5 and 6 depict results that have been obtained using this invention.

FIG. 1 depicts a representation or model 10 of a standard foil, such as an airfoil or hydrofoil, constructed by connecting 200 line segments 11 between adjacent body points 12. This particular airfoil is symmetrical in cross section and has a maximum airfoil section thickness 13 that is 15% of the length of its chord 14. The airfoil chord length has a non-dimensional length of 1.0. Body points are clustered near the leading edge 15 and trailing edge 16 of the airfoil 10 to better resolve the flow at those locations. It will apparent that similar models can also be produced to represent objects with other cross-sectional shapes.

In accordance with this invention there is defined, as an initial condition, a vorticity segment at each body point 12. FIG. 2 discloses one such vorticity segment 20 at the surface of the airfoil 10 centered at location (x.sub.p,y.sub.p). This vorticity segment 20 and others at the surface are infinitely thin. FIG. 2 discloses another vorticity segment 21 off the surface of the airfoil 10 at a location (x.sub.el,y.sub.el). This vorticity segment 21 is prescribed to have the appearance of a flat panel 22 with a finite thickness. Each such segment on the surface of the airfoil 10 carries two velocity generators, namely: (1) a surface vortex distribution lying in the plane of the segment and (2) a potential source. It is assumed that both distributions are uniform over an individual segment. The vortex strength parameter characterizing any panel is the velocity jump or change across the panel. As known, the velocity due to a potential source, .alpha., and .gamma. on a contour C is ##EQU1## where B.sub.n is the expression for the far-field influence of all vorticity segments in terms of a set of coefficients defining the velocity jump across panels as represented by a series of coefficients. The solution of equation (2) is well known in the art.

As will also be apparent, if an object, such as the airfoil 10 in FIG. 1, is defined by N panels, there are 2N unknowns. However, it is also required that the total velocity be zero at the control point, or centroid, of each segment. This constraint produces an equivalent number of equations. The integral of vorticity over a bounded volume is zero when the velocity goes to zero at the bounding surface. Thus the integral of the vortex strength over the body surface must also be zero. The integral of the surface potential source over the body surface must also be zero by continuity. Consequently, there exists a set of 2N+2 equations that enable a matrix solution for the surface quantities through Lagrange multipliers to that the integral constraints are met exactly and the 2N velocity boundary conditions are satisfied in a least-squares sense.

The transfer of vorticity from the surface into the flow is accomplished by the creation of rectangular vorticity elements lying just above the body surface, such a the element 21 shown in FIG. 2. After each element is produced at the surface, it moves to a position directly above the surface and has an initial thickness equal to ##EQU2## and the vorticity of the element is:

where .DELTA.t is the size of the time step, .omega. is the element vorticity, dA is the element area, .gamma. is the surface velocity jump on the panel and dl is the surface panel length.

Each elevated vorticity segment has the same length as the respective segment over which it lies. Therefore, each element is assigned a vortex strength based upon the velocity jump of the underlying surface panel. Its associated velocity field is determined by the Biot-Savart integral: ##EQU3##

Immediately after vorticity is shed in this fashion, minimal surface vorticity is required to satisfy the no-slip velocity boundary condition. As the elements move away from their original position through advection and diffusion, increase surface vorticity, .gamma., is required to meet these boundary conditions until eventually the new vorticity is shed into the flow by creation of a new family of elements. Vortex elements are shed every 0.1 non-dimensional time (t) units.

The evolution of vorticity is prescribed by the vorticity equation. In two dimensions, this equation is: ##EQU4##

The terms on the right side of equation (6) describe the change in vorticity at a point through advection and diffusion, respectively. The effect of advection is accounted for by moving element control points with the local velocity. As shown in FIG. 2, each end point is advected separately allowing for lengthening or shortening as well as rotation of the element. Since the lengths of the elements continually change, the total circulation, defined as .omega. .multidot.dA must remain constant. To achieve this condition, the thickness of the element varies to keep dA constant. Since each vortex segment can infinitely stretch or compress, there is an upper limit on their maximum length. If the vortex segment exceeds this threshold length, it is split into two elements of equal length. Computational this generates a huge number of elements. As an offset, if two vortex segments cross, they are amalgamated in such a way that both linear and angular momentum is conserved. There is also a minimum thickness bound applied on the Kolmogrov length scale that is approximated as: ##EQU5## where .nu. is the non-dimensional kinematic viscosity and .vertline..omega..vertline. is the magnitude of the vorticity in the segment. After this minimum thickness is reached, the total area of the element is allowed to increase with a concomitant decrease in vorticity to satisfy conservation of total vorticity.

The effects of diffusion can be incorporated by the use of conventional random walk techniques that provide a standard deviation of: ##EQU6## where .DELTA.t is the time step size. It is this expression, the placement of the elements as they are shed from the surface, and minimum element thickness that are Reynolds number dependent. In non-dimensional terms, the kinematic viscosity, .nu., is the inverse of the Reynolds number since both chord and freestream velocity are equal to unity.

The stepping in time of the strengths of each vorticity segment and of the segment control points is accomplished using a standard predictor-corrector scheme with one correction applied to the predictor step. Surface pressure was computed according to the method described in Uhlman, J. S., "an integral equation formulation of the equations of motion of an incompressible fluid," Naval Undersea Warfare Center Technical Report 10,086, 15 Jul. 1992. According to that method, that is based on stagnation enthalpy: ##EQU7## The first term in the first interval on the right hand side of Equation (9) accounts for any pitching motion of the air foil 10 or other object. On the surface,

so

When the foregoing steps have been completed, the vorticity of every vorticity segment in an area of interest has been determined. However, as is known, the velocity vector associated with any vorticity segment is dependent not only upon the vorticity of that segment but, to various degrees the vorticity strengths of all surrounding vorticity segments. Conventional processing would involve for the determination of the velocity of any particular vorticity segment, the summation of the influence of the velocity of that vorticity segment as influenced by each of the surrounding vorticity segments. Thus, if there are N vorticity segments, a conventional conversion of vorticity to velocity would include N.sup.2 calculations. This puts a tremendous burden on processing particularly as the resolution with the concomitant increase in the number of vorticity segments required to achieve that resolution.

However, it is also known that the influence of a second vorticity element adjacent a given vorticity segment is greater than the influence of the vorticity of another vorticity segment that is located remotely from the given vorticity segment. Stated differently, the far-field effect of any collection of vorticity segments can be expressed as a multi-pole expansion about the center of a collection of those vorticity segments whose coefficients are merely the sums of the moments of the vortex strengths about that center. This fact allows the construction of a far-field expansion for the collection solely from the knowledge of the vortex strengths and locations. As the number of vortex segments increases, this allows the far-field influence of any collection of vortex segments to be computed to any desired accuracy from a truncation of a multipole expansion.

In accordance with this method and with reference to FIG. 3A, an area of interest concerning the flow over an air foil 10 is defined by a Level 0 square box 20 that is the first level of a tree structure to distinguish vortex segments whose influence must be computed directly from those whose influence may be computed from a multipole expansion. The area of interest bounded by the Level 0 box 20 is defined as containing the complete collection of vortex segments.

The box 20 is then subdivided into a set of four equally sized Level 1 boxes as shown in FIG. 3B that are identified as BOX 2 through BOX 5. Each of these boxes is further subdivided into sixteen Level 2 boxes 22 as shown in FIG. 3C and further divided an array 23 of sixty-four Level 3 boxes as shown in FIG. 3D. If additional resolution is required, additional levels can be produced by further subdividing each successive level. As will be apparent, the number of boxes at any level is 4.sup.L where L is the level from 0 to L.sub.max.

In accordance with this method the coefficients are computed for each of the boxes, BOX 22 through BOX 85 in the Level 3 array 23 shown in FIG. 3D. To obtain coefficients about the center of the box for a truncated multipole expansion of the vortex segments that reside in this box. Once this is accomplished, the coefficients for the truncated multipole expansion for the vortex segments in a box at the next highest level are computed about the center of this box. For example, BOX 9 in FIG. 3C corresponds to BOXES 28, 29, 36 and 37 in FIG. 3D and the coefficients generated with respect to those boxes in FIG. 3D are then used to provide a set of coefficients for BOX 9 in FIG. 3C. This step can be accomplished by simple summation and translation operations. When this procedure has been completed, truncated multipole expansions are obtained for all boxes at all levels.

Once this procedure is completed, velocities for each vortex segment are determined beginning at the second level below the largest box to define an interaction list for each box. This list identifies the boxes for which the interaction must be computed directly. For each of these boxes the coefficients of a Taylor series expansion are obtained about the center and are computed from the coefficients of the multipole expansions and are taken with respect to the center of a smaller box using recursion relationships. The process continues until the Taylor series expansions have been computed for all boxes at all levels.

Once the computation of the Taylor series coefficients is complete, the computation for the interactions begins. This proceeds by computing the direct interaction for all vortex segments in the same box as a field point and for all vortex segments in the boxes on the interaction list of that box. The interaction for all other vortex segments are then computed for the Taylor series expansion for the local box. This procedure is repeated until all desired points and the computation is then complete.

FIG. 4 depicts, in schematic block form, one embodiment of apparatus that can perform the foregoing functions. The apparatus 40 includes a series of input devices, a processing system and an output or display device. More specifically, apparatus 40 includes an object model input 41 that provides a representation of the air foil 10, a flow parameter input 42 and a resolution input 43. The flow parameter input provides information such as free stream velocity, pressure, lift and drag coefficients, Reynolds numbers, kiinematic viscosity and related parameters. The resolution input 43 establishes the maximum value of L and therefore another resolution that will used in providing the velocity representation. A processing system 44 is depicted in this particular embodiment as comprising a series of independent processors. A vorticity segment processor 45 provides the vorticity value for each vorticity segment in the area of interest shown in FIG. 3A. The box overlay processor 46 responds to information from the resolution input for providing the appropriate number of levels of boxes as shown in FIGS. 3A through 3D. A direct velocity processor 47 provides a direct velocity calculation for an individual segment and all the segments in its box and all neighboring boxes while a far field velocity processor 48 provides the approximations of velocity influence from more remote boxes as they represent the sum of all the vortex segments within a box.

A summing processor 50 combines the outputs from the direct velocity processor 47 and far field processor 58 to provide an input to a display processor 51 that can produce any of a variety of outputs on a display device 52.

It will be apparent that this structure can be implemented as a group of discrete processors performing these any other individual functions or by a general purpose computer system utilizing a variety of programs that perform the functions of the various processors shown in FIG. 4.

More specifically, the initial process begins by solving a two-dimensional vortex interaction whereby providing a two-dimensional vortex interaction solution wherein a stream function at a point (x,y) due to point vortex at (.xi., .eta.) is given by

where .GAMMA. is the circulation about the vortex. The velocity components associated with this stream function are then ##EQU8##

The stream function due to a collection of these vortices is given by: ##EQU9##

Although Equations (14) and (15) could be employed to obtain expressions for the velocity components of the system shown in FIGS. 3A through 3D, in accordance with this invention an expansion of the stream function is provided by Equation (16), which yields the following far-field form: ##EQU10## Thus, Equation (16) becomes ##EQU11## It will now be apparent that the coefficients depend solely upon the circulation and location of the point vortices under consideration.

As previously indicated, it is also necessary to move the center of the expansion given above. This can be accomplished by performing a far-field expansion about the translated center (.xi..sub.0,.eta..sub.0) as ##EQU12## When the expansion ##EQU13## is applied, Equation (22) becomes ##EQU14## since s=n+p and t=m+q must remain constant. Thus ##EQU15##

In order to obtain the expansion of the field domain, or velocity, the far-field expansion about the translated center (x.sub.o, y.sub.0) becomes ##EQU16## Using an expansion similar to that used in Equation (22) becomes ##EQU17## and the Taylor series expansion is then given by ##EQU18##

The summation limit L is equal to the maximum value of N less n (i.e., L=N-n). Similar methods can be used to move the center of the expansion given above. Equation (30) then is the expression for the far-field influence of all vorticity segments that is also described with reference to Equation (2).

FIG. 5 plots velocity vectors taken at vorticity element centroids as predicted in accordance with this invention. The sequence shown in FIG. 5 especially shows the adaptive nature of the boxes shown in FIGS. 3A through 3D. Higher concentrations of elements may be found in regions of high vorticity and coherent vortex structures. No vortex segments are indicative of regions of potential flow. This plot is taken for a reduced frequency of .pi./4 taken at different times relative to an oscillation period (t/T). At t/T=0.5, a number of vortex segments may be found near the leading edge. By t/T=0.7, however, the velocity vectors accumulate forming a dynamic stall vortex that may be seen at approximately the mid chord. By t/T=0.9 the vortex reaches the trailing edge and at t/T=0.1, these vortex segments have been shed into the wake. Accumulations of opposite vorticity may be seen nearing the trailing edge forming the trailing edge vortex. These elements advect down stream and eventually form a sinusoidal pattern in the air foil wake.

FIG. 6 shows the stream line plots for a k value of .pi./4 at t/T=0.5 initial formation of a leading edge vortex may be seen the quarter chord region. Unlike experimental results where only one vortex may be seen, this process predicts the formation of two leading edge vortices. The thickness of the separated region looks similar for both the predicted and experimental cases. Downstream there appears a relatively thick separated region consisting of a clockwise vorticity. By t/T=0.7, the vorticity in the separated region is shed into the wake and the leading edge vortices have combined to form a single dynamic stalled vortex which can be seen clearly at approximately 70% of the chord. Just up stream there appear to be smaller vortices in the separated region near the trailing edge. By t/T=0.9 the dynamic stalled vortex is just shed into the wake and is qualitatively similar to experimental results. By t/T=0.1, the dynamic stalled vortex is completely shed into the wake and a counter rotating trailing edge vortex can be seen. In addition the curvature of the streamline suggest an overall counterclockwise circulation about the air foil 10. This has also been observed experimentally. The flow field begins to reattach near the leading edge and continues as t/T=0.3.

These predictions track experimental results. They are obtained, however, without the need for the construction of a model and the physical testing of that model. Changing the shape of an object merely requires changing the model provided to the system. Thus this invention presents an opportunity for obtaining predictions of flow and turbulence around an air foil in an expeditious and cost-efficient manner.

This invention has been disclosed in terms of certain embodiments. It will be apparent that many modifications can be made to the disclosed apparatus without departing from the invention. Therefore, it is the intent of the appended claims to cover all such variations and modifications as come within the true spirit and scope of this invention.