Stream 3D: Computer graphics program for streamline visualization

STREAM 3D: COMPUTER GRAPHICS PROGRAM FOR STREAMLINE VISUALIZATION PETER ELIASSON FFA, Aeronautical Research Institute of Sweden, 161 11 Bromrna, Sweden,

JESPER OPPELSTRUP KTH, Royal Institute of Technology, 100 44 Stockholm, Sweden

ARTHUR RIZZI FFA, Aeronautical Research Institute of Sweden, 161 11 Brornrna, Sweden, and KTH, Royal Institute of Technology, 100 44 Stockholm, Sweden

The present streamline visualization program is based on the technique of mapping the physical domain to a rectangular block. The streamline computation proceeds by integrating the mapped velocity vectors in isoparametric cells and plotting the coordinates of the streamlines thus obtained. The mapping and interpolation use muitilinear functions, and the streamline equations are integrated by the classical Runge-Kutta method with a heuristic strategy for stepsize selection. An example integrating the streamlines in the solution of a vortical flowfieid over a delta wing at angle of attack demonstrates the use of the program.

INTRODUCTION In the study of field theory physicists often use the concept of field lines to improve their understanding. Michael Faraday introduced the concept in the early nineteenth century as an aid in visualizing electric, magnetic, and gravitational fields. A line in such a vector field is an imaginary curve drawn so that its direction at any point is the same as the direction of the field at that point. Since at any one point the field can have but one direction, only one line can pass through each point of the field, and thus the curves never intersect. Today, in addition to the study of electricity and magnetism, the concept is also very useful for visualizing velocity, vorticity, or surface stress fields of fluid flow, particularly those fields determined by numerical methods at discrete points in space. If the vector field is a force field, the field lines are known as the lines of force, and Accepted May 1989. Discussion closes April 1990.

162

Adv. Eng. Software, 1989, Vol. ll, No. 4

if the field is the velocity field of a flowing fluid, they are known as streamlines. Although the concept applies generally to any field, we develop it here in the context of flow fields that have been determined numerically. Consider the point P of the vector field f in Fig. 1. The field line through the point is a general space curve. The position vector r gives the location of P as a function of some parameter t that varies along the field line. It may be time, or arc length, or any other convenient measure. The tangent to the curve at P then is dr/dt and thus from the definition of the field line we derive the equation dr - - = f(r) dt

(1)

which determines the line. By solving this equation at a given instant of time we can construct a picture of the field lines. If the field depends on time the picture will change from instant to instant, but if it is steady, the picture remains the same for all times. When the velocity field of a flowing fluid is steady, the streamlines then become identical to the so-called path lines which are the curves traced out in space by the moving fluid particles*. A picture of the streamlines helps us to see the flowfield and therefore plays an important part in the analysis and understanding of fluid motion. If the vector field f is continuous, the integration of Equation (1) proceeds straightforwardly without additional consideration. Here we are concerned with the integration of Equation (1) when f has been previously determined by some numerical solution at discrete points in space, i.e. the field is discrete, and further consideration is necessary. In the following sections we

* Consult Prandtl and Tietjens ~ for a detailed discussion regarding the relation between path lines, streamlines, and the so-called streak lines.

~ 1989 Computational Mechanics Publications

T[-~i field fine

P

~'(t÷dt)

.,,~I

f

f

?(t)

f

STREAMLINE Figure2. Field line through a computational cell defined by the eight corner points of a discrete vector fieM Figure 1. Geometry underlying the tangency condition dr/dt = f(r) o f the the fieM line

value of f(rp). The interpolation must have finite support over the given discrete values of f, and we must have some strategy to move this support along with the integration as it proceeds down field. A convenient strategy is to define the interpolation over the hexahedral cells of the mesh (see Fig. 2). The discrete values of f are given at the eight vertices (corners) of the cell and one obvious choice is to use these as the support for trilinear interpolation. It is also the natural choice to decide where to locate the support. One need only know which of the six cell faces the streamline crosses in order to determine in which new cell the integration should proceed. Using this strategy the program then carries out the integration of Equation (2) together with the graphical representation of the field lines in a rather straightforward application of the concept of isoparametric trilinear finite elements. Our discussion here treats the method applied to a specific problem in aerodynamics - the tracing of streamlines around an aircraft wing by integration of

describe why the treatment of discrete values of f leads naturally to a transformation of Equation (1) to an isoparametric space. We then present a numerical method to integrate the transformed equation, and conclude with an example of steady streamlines calculated in a vortex flow over a delta wing. FORMULATION

OF THE

DISCRETE

PROBLEM

The field lines are the solutions r ( t ) - [ x(t), y ( t ) , z (t)] T of the ordinary differential system Equation (1) dr dt

-- = f[r(t)]

(2)

but now f may not be known at the field line point rp. So we have to devise some means to interpolate the

Z

J

¢x Figure 3.

f

t>Y

U

X

Mapping L from the ~, ~1, ~ computational space to the x, y, z physical space Adv. Eng. Software, 1989, Vol. 11, No. 4

163

the velocity components of the computed flowfield. The input to the program is a grid and a velocity field computed previously on the points of the grid. These input data are obtained in the following way. A grid generation procedure computes the points of the grid from the shape of the wing. Then the velocities are computed according to the theory of inviscid flow by the program W I N G A 2 (see Ref. 4). The field line computation assumes the vectors are at grid points, i.e. at the vertices of the computational cells. If the velocity vectors are computed at the cell centers, as is the case with the centred finite volume method of W I N G A 2 , a preliminary interpolation must be performed to obtain the vectors at the cell corners. The mesh in our example is organized as a large three dimensional array of N I × N J × N K points which are corners of ( N I - 1) x (N J - 1) × ( N K - 1) hexahedral ceils. The mesh generation implies a mapping from the computational (4,~/,~')-space to the physical Ix, y, z] -space. For field line computation we only need to use this mapping locally in each individual cell (see Fig. 3), i.e. r=

= L(p); where p =

Thus L can be approximated from the coordinates of the cell corners, and, in this way, any mesh generation procedure can be used without changing the streamline program.

SOLVING THE DISCRETE PROBLEM

Isoparametric mapping technique To solve the discrete system (2) we must define the field at any point of space from the grid data. Clearly the mesh transformation is one interpolation scheme, but one which is too cumbersome to use in its nonlocal form. The simplest interpolation procedure is a trilinear function over each cell, and a cell-by-cell procedure immediately comes to mind. That is, the streamline is computed by integration of (2) inside a cell, and when the line reaches the boundary of the cell, it continues into the neighbor, where a new interpolant for the field vectors is computed and the calculation then proceeds. Depending on how the grid mapping L is determined, a trilinear interpolation can be identical to L. The next question is how to handle the definition of the cell. As will be clear from the sequel, using the isoparametric technique is more economical than computations entirely in the physical space, because interpolation over, and checks for containment of, a cubic cell are much simpler than for a general hexahedral cell defined by its vertices. We then have to specify the streamline equation in the mapped space. The field lines in computational space are traces of the solutions to d dt o = g [p(t)]

(3)

where g ( o ) = J - l ( p ) f [ L [ o ( t ) ] ] and J=OL]O0 is the Jacobian of the mesh mapping. There are several conditions that the mapping must satisfy. First of all, we can 164

Adv. Eng. Software, 1989, Vol. 11, No. 4

only calculate a streamline at a point if the field l' is single-valued there. The types of meshes in our application contain singular points of either polar or parabolic character. That means that a cell edge can degenerate to a point, or a cell face can collapse to a line or a point at such a singular point. In other words, there can be several grid points occupying the same location in space. The flow solver, however, ensures that although there may be multiple values of f, they are all identical. The Jacobian formally must be non-singular, even at singular points in the mesh. This depends on how one chooses to approximate OL/Oo. Now we would like to define an interpolation scheme for the right hand side g(p) that has enough continuity to guarantee an (almost) classical solution p(t) to Equation (3). If we consider only consistent approximations in the sense that J = OL/Op is exact, continuity of g is too much to ask for, since that would require L to be in C x, which is practically impossible because of the computational cost. As shown below, using another (inconsistent) approximation to J allows construction of a continuous g within the cell. However, it is possible then that g suffers a discontinuity only at cell faces; and so we have to insure that a line crossing the face into one of the cells continues into that cell for some small distance, i.e., sign(g, n ) l = sign(g, n) N

(4)

where the superscripts refer to the two neighboring cells sharing a face with normal n. It turns out that the use of a piecewise (cell-by-cell) trilinear L satisfies this condition if the cells in the physical space are not too distorted. It is clear from the construction by interpolation that L will be continuovs, but will have discontinuous derivatives.

Algorithm to compute a field line With the above considerations the algorithm can now be summarized as follows. Let rijk, i = 1..... 741, j = 1..... N J, k = 1..... NK, be the cell corners with associated field vectors fijk, and let Cijk be the cells numbered such that the corner closest to the origin of C~ik is r~yk. The computational coordinate system (~,~, ~') of a cell is so chosen so that 4, corresponds to the/-direction, ~/to the j-direction, and ~"to the k-direction. Then: a) Let the field line start at (40, ~0, ~'0) in Cijk. b) C o m p u t e the coefficients of the mapping L,jk c) Compute the J, and the mapped field vectors qs = J s - l f s at the eight cell corners. d) Integrate the field equations using the trilinear interpolant to the gs until the field line reaches a face of C~j.k, say the one with ~- = ~'* = 1, with 4 = 4", r/= ~/ . Output the physical coordinates Lijk(4*, ~/*, ~-*) for plotting. e) If the face is a boundary face, i.e. k = N K - 1, stop, else, set k = k + 1, 40 = 4", ~/0 = , * , ~-0 = O, and go to step a) and repeat. In steps d) and e) the obvious modifications are made for the five other cases when the line reaches some other face of the cell. We now describe the major parts of the algorithm.

Constructing the Mapping Lug and Jacobian J The mapping Lijk itself is used for the computation of the computed points on the streamline in physical space. It is also needed as the interpolant of the field vectors f~jk which together with J yields the field vectors g in the computational space. The eight corners of the computational cell Cug are numbered i + I, j + J, k + K where I, J, K = O, 1 and the computational coordinates of corner IJK are, of course, (~, r/, ~ )IJK = (I, J, K )

Defining the one-dimensional blending functions

because J --~ 0L/0o. Another approximation Jel = (ri+ l,j,k - ri, j,k ) Je2 = (ri, j+ l,k -- ri, j,k )

Je3

=

(ri, j , k + l - ri, j,k)

does not have this continuity property, but the result is exact (i.e. J = 0L/0p) since L is a trilinear interpolant, hence also J is a linear function of distance (i.e. a coordinate) along an edge where the other two coordinates are constant. We have used the second approximation in the examples below, but the first also gives comparable results.

~b0(t) = 1 - t; k l ( t ) = t we have that ~bz(J) = 61,j for I, J = 0, 1. The tensor product basis functions

~ u ~ = ~bt [ ~] ~ j[~ ] ~ [ ~'] then satisfy the interpolation relations ~)IJK(L, M , N ) = ~IL~JM¢~KN

Then the trilinear interpolant Lij,($, ~, ~) to the corner data hlJl( = {rg} IJK becomes Lijk [h(~, ~/, ~')] = C t j r h l j ~ l j ~ ( ~ , ~, ~) which may be computed in 3 x (8 x 3 + 7) = 93 floating point operations if expanded as is. But if L is to be evaluated m a n y times, as is the case here, it pays to preprocess the expression for L to the form L = c + a l ~ + a 2 r / + a3~" + bl~'~/+ b2~~ + b3~7 + d~B~"

Integration step size Once the interpolation of the velocity vector is done, the integration can proceed. The integration in a cell is done by the fourth order R u n g e - K u t t a method and carried out in the ~, ~/, ~- coordinate system from the first value in the new cell to where the integrated streamline leaves the cell. The step length is calculated for each cell so that approximately N steps are taken in the cell where N is set by the user. If we determine an average velocity by 8 Uave = 1/8 Z [ g ' l i=1

then a step size At is given by At = l/(vave N ) . Eliasson 2 has experimented with different numbers of steps N, and he found convergence to full plotting accuracy of the streamline integration for vortex flow around a delta wing with a value N = 3.

evaluated as, say L = c + ~'(a3 + b l ~ 4- ~(b2 4- dr)) + ~(al

+

b3r/) + az~

by computing c, ai, b i and d from hts~. In this case the work becomes 42 flops not including the preprocessing. We get c = hooo

;hi = hoxl - c -

a2-

Program structure The algorithm is carried out by calls to a sequence of sub-routines as indicated in Fig. 4. Table 1 lists the functions of these routines.

> COEFF

a3

> JACOBI

a~ = h~oo- c;b2 = h101 - c - ax - a3 a2 = h o ~ o - c;b3 = h l ~ o - c a~=hoo~-c;d

>

> SOLVE

a~ - az

>F

=hlll-c-ax-a2-a3-bl-b2-b3

Notice that the same procedure is used to compute the field vectors gt in computational space and the points on the field rt lines in physical space. The hlj~ in the latter case are the physical coordinates ri+1,j+J,k+r of the cell corners. The Jacobian J is needed to obtain g. The Jacobian at a corner of a cell is computed by differences over the adjoining edges of the cell. At the corner rijk, since the cell sides are of unit length in the computational space, we can obtain the three columns of J Jel

> USTART

> COEFF

> UVEL

>

MAIN

>

>l

> COEFF

JACOBI

> STEPL

= ( r i + l,j,k -- r i - 1 , j , k ) / 2

>F

Je2 = (r/,j+ l,k -- Fi, j - i , k ) / 2 Je3 = (ri, j,k+ l - ri, j . k - t)/2 where ei are the unit base vectors. This is an approximation that results in both L and J being continuous across cell faces, but it is not exact

> SOLVE

> INTEG > CUT

> NEXT

Figure 4.

Calling sequence o f the routine

Adv. Eng. Software, 1989, Vol. l l , No. 4

165

Table 1. The subroutines o f the algorithm Subroutine:

COEFF: SOLVE: JACOBh UVEL: INTEG:

F:

CUT: NEXT: STEPL:

USTART:

Function:

Calculates the 24 coefficients for the transformation r = L(p) and the interpolant g(o) Solves a third degree equation system. SOLVE is used to determine the eight velocity vectors in the ( ~ ) coordinate system Calculates the Jacobian J Calculates the eight velocity vectors g and calls for COEFF to calculate the interpolant g(o) Integrates in a cell from po to pX using the fourth order Runge-Kutta method. When the integrated streamline has advanced to outside of the cell, INTEG calls for subroutine CUT to cut the streamline at the cell edge Calculates r = L(p) or g(p) L[g(p)l Cuts a line at the cell edge Determines the start coordinates for the next cell and the cell indices 1, J, K for that cell Calculates the step length for one cell so that approximately N steps will be taken in the cell. N is set by the user Calculates the start value in the ( ~ ' ) coordinate system and the cell indices from the given start value r in the (xyz) coordinate system. To calculate p from r = L(O) a nonlinear system of equations are solved using Newton Rapbson's method.

F o r the f o u r t h - o r d e r R u n g e - K u t t a m e t h o d , f o u r e v a l u a t i o n s o f the i n t e r p o l a n t g ( p ) are m a d e per step. O n e e v a l u a t i o n o f g ( o ) costs a b o u t 60 o p e r a t i o n s . T h e f o u r e v a l u a t i o n s plus s o m e e x t r a o p e r a t i o n s m a k e a b o u t 280 o p e r a t i o n s per step. If three steps are m a d e in a cell, the i n t e g r a t i o n requires a b o u t 850 o p e r a t i o n s per cell. T h a t is a p p r o x i m a t e l y the s a m e w o r k t h a t is r e q u i r e d for the i n t e r p o l a t i o n o f the velocity vector. The t o t a l n u m b e r o f o p e r a t i o n s per cell is t h e r e f o r e a p p r o x i m a t e l y 1650 per cell if an a v e r a g e o f three steps per cell are taken.

E x a m p l e o f vortical streamlines The s i m u l a t i o n o f vortex flow shed f r o m the edge o f a delta wing has been the o b j e c t i v e o f a recent international s t u d y 3. The wing in that s t u d y is a delta with a 65 deg. leading edge sweep a n d a c r o p p e d side edge. W e have p r o d u c e d solutions to b o t h the Euler a n d N a v i e r Stokes e q u a t i o n s for flow M = = 0.85 a r o u n d this wing at c, = 20 deg. angle o f a t t a c k . Figure 5 presents two views o f the streamlines i n t e g r a t e d f r o m the Euler solution with S T R E A M 3 D . Because the vortex is shed a l o n g the leading edge, we begin the s t r e a m l i n e i n t e g r a t i o n at a n u m b e r o f positions a l o n g the leading edge. The oblique view (Fig. 5a) and the side view (Fig. 5b) clearly show the vortical m o t i o n o f the streamline. T h e velocity field here has been c o m p u t e d with a mesh o f a b o u t 300,000 grid points, a n d the i n t e g r a t i o n o f these streamlines t h r o u g h the field was carried out in less t h a n 1 1/2 rain o f C P U time on the C Y B E R 205 supercomputer.

FIELD LINES IN A SURFACE A c o m m o n w a y o f visualizing fluid flow in the l a b o r a t o r y is by s m e a r i n g d y e d oil o n t o the solid surfaces e x p o s e d to the flow. T h e streaks in the p a t t e r n which d e v e l o p indicate the d i r e c t i o n o f m o t i o n o f the fluid in the b o u n d a r y layer. T h e r e are two c o u n t e r p a r t s to this p h e n o m e n a in n u m e r i c a l l y s i m u l a t e d flowfields. In an Euler s o l u t i o n it is the streamlines that lie in the surface o f the wing, because the b o u n d a r y c o n d i t i o n s ensure that the velocity is parallel to the wing. In a N a v i e r - S t o k e s s o l u t i o n the b o u n d a r y c o n d i t i o n s require that the velocity vanish on the surface. N o w it is the viscous forces t h a t d o m i n a t e . O n e resorts then to the n o t i o n o f the limiting streamline, i.e. the streamline infinitesimally a b o v e the surface (see Refs 6 a n d 7). A limit process is involved, but it turns out that the

a) oblique view

b) Figure 5.

166

side

view

Streamlines c o m p u t e d f r o m an Euler solution M= = 0.85 ~ = 20 °

A d v . Eng. Software, 1989, Vol. 11, No. 4

i

Figure 6. Skin-friction lines c o m p u t e d on the upper surface o f the delta wing. Laminar N a v i e r - S t o k e s solution M~ --- 0.85 ~ = 20 ° Re = 2.38 x 106

limiting streamline is equivalent to the so-called skinfriction line. In this case the field is the vector component of the stress tensor r that lies in the surface, the socalled skin friction force r~. In either case the resulting surface field lines can be integrated by a twodimensional variant of the algorithm we have presented. If the surface has the parameter representation r = S(O), where now O = [,~] the in-surface component o f a vector is the orthogonal projection onto the range of the Jacobian J = OS/Oo, then the vector in parameter space g=

v'

may be computed as g = [ j T j ] - I j T f .

The same technique can also be used to draw skin friction lines by integrating the vector r • n (the surface normal vector component of the stress tensor r) instead of the projected velocity vector.

Mapping

The method to compute surface streamlines or skinfriction lines is closely connected to the generation of space streamlines. A multilinear transformation is set up between the three dimensional solid surface of each computational cell and a two dimensional coordinate system where interpolation and integration are performed. The mapping between three dimensional physical (x, y, z)-space and the two dimensional computational (~,~)-space is set up as before. The transformation r = S(O) now is written in the form r = S ( p ) --- al + a2~ + a3~7 + a4~7

The twenty-four coefficients of at, a2, a3 and a4 are determined in corresponding fashion as for space streamlines. After the computation of the Jacobian and the projected vectors g of each corner of the cell face, the interpolation and integration can proceed as for space streamlines. The bilinear ansatz for the vector g(~, ~7) in the two-dimensional space of the surface is g(~, ~7) = bl + b2~ + b3~ + b4~7

E x a m p l e o f skin-friction lines

Here we begin with our Navier-Stokes solution ~

computed for flow Moo = 0.85 around the 65 deg. delta wing at 20 deg. angle of attack and Reynolds number Rec = 2.38 x 10 6 based on the root chord. Figure 5 presents the skin-friction lines integrated by our technique. The integration begins at points distributed along chord-wise lines from the leading to the trailing edge. Converging skin-friction lines signify a region of flow attachment, and diverging lines signify flow separation. The diverging lines in Fig. 6 indicate the secondary flow separation under the primary vortex 5.

CONCLUDING

REMARKS

We have presented a technique to integrate the field lines of a vector field that is given by a numerical simulation. The integration is done by a fourth-order R u n g e - K u t t a scheme in computational space. It makes use of an isoparametric mapping from the real space to the computational space because the integration is more efficient there. Interpolating in physical space requires only the computation of the 24 coefficients in the interpolation dr/dt = L [f(r)]. But since the interpolation is not made over the simple unit cube, this leads to three systems of equations of degree eight for the coefficients, instead of the eight systems of third degree as we have. Determining where the streamline crosses one of the six faces of the cell in physical space is also more complicated. The method is presented for the general threedimensional case, and an example of streamlines around a delta wing demonstrates the use of the computer program. The same technique can be specialized to the problem of field lines in a surface. An example of skinfriction lines computed f r o m a N a v i e r - S t o k e s solution demonstrates this variant of the program.

REFERENCES

1 Prandtl, L., and Tietjens, O. G. Fundamentals of Hydro- and Aeromechanics, Dover, New York, 1967 2 Eliasson, P. Graphical Presentation of Streamlines in ThreeDimensions, FFA TN 1986-53, Stockholm, 1986 3 Elsenaar,A., and Erikson, G. (eds), Proc. Intl. Vortex Flow Experiment on Euler Code Validation, FFA Rep. Stockholm, 1986 4 Rizzi, A., Drougge, G., and Purcell, C. J. Euler Simulation of Shed Vortex Flows over the 65 Degree Delta Wings, in Ref. 3, pp. 289-343

A d v . Eng. Software, 1989, Vol. 11, No. 4

167

Rizzi, A., Drougge, G., and Moiler, B. Navier-Stokes and Euler Solutions for Vector Flow over a Delta Wing, in Proc. Symposium Transonicurn 111, Springer, 1988

168

A d v . Eng. S o f t w a r e , 1989, 17ol. 11, No. 4

6 7

Perry, A. E., and Fairlie, B. D. Critical Points in I'[ov, Patlc)ii,, Adv. Geoph_vsics, 1974, 18, 299 315 Dallmann. U.: Topological Structures of Thrcc-l)mlensional Flow Separation, AIA.4 Paper 1983, 83-1735