Below are the first 10 and last 10 pages of uncorrected machine-read text (when available) of this chapter, followed by the top 30 algorithmically extracted key phrases from the chapter as a whole.

Intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text on the opening pages of each chapter.
Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.

Do not use for reproduction, copying, pasting, or reading; exclusively for search engines.

OCR for page 683

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
SESSION 13
LIFTING-SURFACE FLOW: UNSTEADY VISCOUS METHODS

OCR for page 683

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
This page in the original is blank.

OCR for page 683

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Direct Numerical Simulation of Unsteady Propulsor Blade Section Forces on a Parallel Processor
W.B.Coney,1 S.R.Breit,2 and J.R.Webb1
(1BBN Systems and Technologies, USA, 2Kendall Square Research, USA)
ABSTRACT
This paper describes a numerical method, which has been implemented on a parallel processor, for the direct simulation of inviscid unsteady blade section forces, and its application to several two-dimensional blade-gust interaction problems.
The parallel implementation is discussed in some detail. Timing results obtained on a parallel computer are presented and nearly linear parallel speed-ups are demonstrated.
The application of the present method to three model problems is discussed and results are compared to theoretical, experimental and numerical solutions.
INTRODUCTION
The time-harmonic (blade-rate and its harmonics) and stochastic unsteady forces produced by interaction with non-uniform inflow are the dominant components of the low-frequency acoustic signature of a non-cavitating marine propulsor. Much research has been devoted to analytical and approximate numerical methods for predicting these forces and there has been considerable success with these approaches. However, with the recent dramatic increases in the speed of computers and the advent of parallel computers, direct numerical simulations of unsteady flows through marine propulsors are becoming increasing feasible. The ultimate benefit of such direct simulations will be the ability to deal with arbitrarily complex propulsor geometries and to include possible interaction between the blade-rate and turbulence ingestion mechanisms which are presently treated separately.
The prediction of the response of a two-dimensional blade section to an unsteady gust can be regarded as a preliminary step toward the fully three-dimensional model of the propulsor subject to an unsteady inflow. The force developed in response to harmonic gusts of varying nature has further been the subject of much analysis, leading to both analytical, Sears [1], Naumann and Yeh [2], Atassi [3], etc., and linearized numerical, Atassi and Akai [4], Chiang and Fleeter [5], Hsin [6], etc., approaches. Unsteady blade-gust interaction problems therefore provide a natural test bed for the exploratory development of a numerical simulation of unsteady propulsor forces.
This paper describes a numerical method, which has been implemented on a parallel processor, for the direct simulation of inviscid unsteady blade section forces, and its application to several two-dimensional blade-gust interaction problems. The results selected for presentation highlight circumstances for which “classical” unsteady airfoil theories, such as that of Sears [1], and numerical procedures based on a linearizing assumption of potential flow, for example those of Hsin [6] and of Chiang and Fleeter [5], appear to adequately predict the unsteady force, and circumstances for which they cannot. Three model problems are investigated:
transverse gust interacting with a foil
pitching foils upstream of a stationery foil
distributed vortex interacting with a foil.
The first problem is well modeled by an analytic solution, the second could at least potentially by modeled by an unsteady potential flow simulation, the third, however, requires use of a method which preserves detail of the vortical flow.

OCR for page 683

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
NUMERICAL APPROACH
In choosing a numerical method for direct simulations, it should be recognized that the unsteady flow through a typical marine propulsor is characterized by:
very low Mach number (below 0.01),
high Reynolds number (above 106),
complex geometry, and
small unsteady pressure (compared to the mean pressure),
or small unsteady thrust (compared to the mean thrust).
With these issues in mind, we chose to implement the time-accurate pseudo-compressibility (TAP) method with upwind-biased spatial differencing, as described by Rogers & Kwak [7] for our simulations. In the limit where the pseudo-time step size is infinite, this method is nearly the same as the Newton iteration scheme of Athavale & Merkle [8]. The TAP method provides an incompressible flow solution which is a very good approximation of a low Mach number flow. In contrast, methods which account for compressible effects break down at very low Mach number. Our initial experience in applying this method to a two-dimensional model problem, with focus upon one aspect of the blade-vortex interaction, is described in Breit, et. al. [9].
The time integration scheme is second-order accurate and fully implicit, including the boundary conditions. Viscous terms are, for the present, neglected. A sparse, block-tridiagonal, linear system must be solved at each pseudo-time iteration. Unlike Rogers & Kwak [7], who use a line relaxation scheme to solve the linear system, we use the approximate LU factorization scheme recommended by Merkle and Athavale [8].
PARALLEL IMPLEMENTATION
Implementing the approximate LU factorization scheme on a massively parallel computer is a programming challenge because the degree of parallelism varies greatly during the course of the factorization procedure. The parallelization was accomplished on a BBN TC2000™ parallel computer. The implementation relied to a large extent on the high complete-exchange bandwidth and high random-access bandwidth of the TC2000 computer.
In theory, the highest parallel efficiency can be achieved if the smallest computational task that can be done in parallel involves the computations for a single grid point. We refer to this as a point parallel approach as opposed to a domain decomposition approach, where the smallest task involves computations for a sub domain of the grid. The domain decomposition approach requires less interprocessor communication than the point parallel approach, but at the cost of reducing the maximum achievable parallel efficiency. The choice of which approach to implement on a particular massively parallel computer depends on the characteristics of its interprocessor communication network. The interprocessor communication is fast enough on the TC2000 to enable a point parallel implementation of the LU algorithm.
The parallelization process was simplified by exploiting the memory management extensions and parallel directives provided by TC2000 FORTRAN. Some data are shared using an interleaving mechanism, with data declarations of the form
common, shared, interleaved /foo/ a(100),..
In this example, the set of 16-byte cache lines that contain array a are dynamically interleaved among the memories of the processors. Thus interprocessor communication is effected through shared memory instead of explicitly passing messages or packets between processors. Automatic load balancing is achieved by dynamically allocating computational tasks to processors. This approach avoids special mapping strategies where tasks and data are preallocated to specific processors. One benefit of this approach is that the number of lines of code that are associated with the parallel implementation is small compared to the total size of the program.
The LU method involves a forward elimination sweep over the grid and then a back substitution sweep in the opposite direction. Consider only the forward elimination sweep and denote a point in a structured 3D grid by the indices (i,j,k). The elimination for a point (i,j,k) can not be done until the points (i-1,j,k), (i,j-1,k) and (i,j,k-1) have been completed. Thus the elimination must start at the point (1,1,1), then proceed to the points (2,1,1), (1,2,1), (1,1,2), and so on. Obviously the degree of parallelism varies widely during the course of the forward elimination sweep and there is a corresponding variation in the degree of parallelism during the back substitution sweep. It is therefore impossible to achieve linear speed-up of this algorithm. The maximum speed-up depends on the problem size and the number of processors.
In order to parallelize the approximate LU factorization we define a set of grid points in some plane n which satisfy
i+j+k=n+2
(1)

OCR for page 683

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
For example, on an N×N×N grid, there are 3N-2 such planes and the number of grid points in each plane is given by
(2)
Such a plane is illustrated in figure 1 for a 16×16×50 grid (used in the investigation [10] of a three dimensional vortex convecting through a converging/diverging channel). In the forward elimination step of the LU algorithm, we loop sequentially over n in ascending order. In our N×N×N example, from 1 to 3N-2, but the operations can be done in parallel on points which lie on the same plane as defined above.
Figure 1. Planes of grid points available for point parallel calculations in the elimination/substitution sweeps of approximate LU factorization algorithm.
The parallelization is effected by means of the TC2000 FORTRAN spread do command, with the do statement acting over the number of grid points, n, on the current plane. The spread do potentially executes each loop iteration on a separate processor.
spread do
do m=1, n
…
end do
end spread
It will be noted from Eq. 2 that for the cubic grid the number of points on the planes range between 1 and 3/4N2, thus the degree of parallelism varies during the course of the forward elimination sweep. The back substitution sweep is similar except that we loop sequentially over n in descending order. Due to the variation in the degree of parallelism during the course of the forward and backward sweeps, the maximum achievable parallel efficiency decreases as the number of processors increases. The present approach is by no means the optimal one, but it does produce good results on a moderate number of processors and was easy to implement
Parallel Timing Results
Figure 2 shows the parallel speed-up obtained with different numbers of processors for an implementation of an approximate LU factorization, similar to that discussed here. These particular results are for the parallelization, performed by a BBN team [11,12], of a parallel benchmark released by NAS. This benchmark is part of a suite of such benchmarks which are intended to, as a whole, “mimic the computation and data movement characteristics of large scale computational fluid dynamic (CFD) applications [13].” The grid size for this calculation was 643.
Figure 2. Parallel speed-up achieved for approximate LU factorization on a BBN TC2000 parallel computer.
The solid line in Fig. 2 represents the speed-up for a computation of ideal parallel efficiency, η,
.
(2)

OCR for page 683

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Note in Fig. 2 that nearly ideal speed-ups are achieved for as many as P=32 processors, with quite good speed-up (η=91%) obtained for the maximum number of processors, P=55, utilized in this study.
The fall-off in parallel efficiency with number of processors can, in part, be attributed to the previously discussed decrease in the maximum achievable efficiency for the present algorithm due to the varying degree of parallelism throughout the approximate factorization process. An additional drop-off in efficiency is due to the use of the relatively large number of synchronization barriers which are needed to achieve the parallelization as presently implemented. It should be noted that both sources of parallel efficiency loss would play a smaller role as the problem size is increased.
Figure 3. Parallel timing results achieved for present implementation on a time-shared, multiuser system. The various symbols denote runs with differing grids or boundary conditions.
Figure 3 shows some parallel timing results obtained with the present code. Results are for various unsteady simulations performed on the 2-D grids of Fig. 4, Fig. 6 and Fig. 7. Elapsed time per computational time step is plotted against the number of processors used for the computation. The Fig. 3 results are particularly noteworthy in that they are obtained for unsteady calculations using real geometries on a time-shared, multi-user computer, simultaneously in use by other users. This is in contrast with the more usual benchmark timing results, typically obtained under ideal operating conditions, on a “clean” computer. The present results reflect practical speed-ups obtained by performing the computations in parallel. The speed-up is noted to be very nearly linear with the number of processors, up to the largest number of processors utilized in this study.
Even though the BBN TC2000 is no longer in production, we believe that these results illustrate the potential of parallel computing for the practical solution of computational fluid dynamics problems.
RESULTS
Foil in Steady Flow
Figure 4 shows a grid used in the present study. The section shown in Fig. 4 is here denoted “Foil 1”. Foil 1 is a NACA 16 section with a thickness/chord ratio of 11.56% and camber to chord of 3.36%, and an a=0.8 mean camber line. The angle of attack was set to 0.8 degrees. The grid shown in Fig. 4 is a C-grid with 281×71 grid lines.
Figure 4. Grid used for Foil 1 transverse gust calculations, a C-grid with 281×71 grid points. For clarity of illustration only every second grid line is shown.

OCR for page 683

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 5 gives the steady-state pressure distribution on the surface of Foil 1 obtained by the present method and with a boundary element method [6]. We note that the differences are small, occurring mostly in the vicinity of the trailing edge. The lift coefficients, CL=L/(1/2ρV2c), obtained with the two methods are nearly identical, CL=0.570 from the present method, vs. CL=0.576 from the boundary element method.
Figure 5. Steady pressure distributions computed for Foil 1 by present method and by a boundary element method [6].
.
Transverse Gust
Transverse gusts interacting with three foil geometries were investigated. These geometries were:
Foil 1: NACA 16 thickness form, a=0.8 mean camber line, t/c=0.1156, f/c=0.0336, α=0.8 degrees.
Foil 2: thickness form identical to that of Foil 1 without camber or angle of attack (see Fig. 6), t/c=0.1156, f/c=0, α=0.
Foil 3: one-half the thickness of Foil 2 (see Fig. 7), t/c=0.0578, f/c=0, α=0.
The grids used in the transverse gust computations are shown in Fig. 4, Fig. 6 and Fig. 7. In each case a C-grid with 281×71 grid lines was utilized. The grids were generated with the aid of the INMESH [14] grid generation software.
Transverse sinusoidal gusts were modeled by imposing vertical velocity perturbations on the inflow boundary, the grid line which defines the outermost “C” in the grids of Fig. 4, Fig. 6 and Fig. 7. The imposed vertical perturbation velocity on the boundary was given by:
(3)
Here k is the reduced frequency, k=ωc/2V. A gust amplitude, , of 4 percent of the steady inflow velocity, V, was employed in the present study. Reduced frequencies of 0.5, 1.0, 2.0 and 5.0 were investigated. A time step size of ΔtV/c=0.025 was used in all of the presented results.
Figure 6. Grid used for Foil 2 transverse gust calculations, a C-grid with 281×71 grid points. For clarity of illustration only every second grid line is shown.

OCR for page 683

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 7. Grid used for Foil 3 transverse gust calculations, a C-grid with 281×71 grid points. For clarity of illustration only every second grid line is shown.
Figure 8 shows the resulting amplitude of the lift response to the transverse gust as a function of the reduced frequency, k. The unsteady lift amplitude in Fig. 8 is normalized by the quantity , the expected quasi-steady value. Also shown in Fig. 8 is the classical solution due to Sears [1] for an airfoil of zero thickness. We note that agreement between the present method and the Sears solution is relatively good. The value of the phase was also checked against the theoretical value for several of the investigated cases, and was similarly found to be in good agreement
Flapping Foil
Detailed measurements of the flowfield about a foil experiencing unsteady time harmonic gusts have recently been performed at the Marine Hydrodynamics Laboratory at MIT [15,16,], and have been the subject of an ONR sponsored workshop at which comparisons between the experimental results and numerical simulations were presented [17]. The gusts were created by placing two small “flapping” foils in the water tunnel upstream of a larger stationery foil about which the measurements were made. This arrangement is shown schematically in Fig. 9
Figure 8. Unsteady lift amplitude due to a transverse gust. Results are shown at several reduced frequencies for the three foil geometries investigated with the present method. The classical solution due to Sears [1] is also shown.
Figure 9. Schematic representation of the flapping foil experiment installed in the MIT water tunnel.
The MIT experiment focused upon measurement of the unsteady velocities in the boundary layer of the stationery foil, with particular emphasis on the trailing edge flows. While the present code is inviscid, and therefore cannot capture the character of such flows, the MIT data still provides a challenging simulation problem for which at least some comparisons with experimental results can be made.

OCR for page 683

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 10. Grid used for flapping foil simulation, a C-grid with 251×51 grid points.
The MIT researchers measured both steady and unsteady velocities at points defining a box surrounding the stationery foil. It was originally our intent to simply extrapolate this data to provide an inflow boundary condition for the numerical simulation. However, there proved to be insufficient data to so proceed. Therefore, a far field boundary condition model was constructed with potential flow singularities to model the effects of the primary elements of the flapping foil experiment. This model is described in the Appendix. Parameters in the model were estimated using a non-linear least squared error approximation to unsteady data and a direct calculation of the steady circulation.
The foil geometry used for the flapping foil simulation is similar to that of the previously described Foil 1, and corresponds to the “unbounded” foil shape which the MIT researchers were attempting to experimentally model in the water tunnel. A somewhat different foil geometry was actually tested by MIT in an attempt to make up for tunnel wall effects. The angle of attack and maximum camber used in the simulation were reduced somewhat so that the predicted lift matched the experimental value. This reduction in lift is attributed to viscous effects which are not modeled by the present method. Fig. 10 shows the grid used for the numerical simulation. We note that the inflow boundary of this grid is quite close to the leading edge of the stationery foil, in order to keep the flapping foils outside of the computational domain.
Figure 11. Predicted contours of vorticity induced by the flapping foils.
Figure 12. Unsteady lift, , induced on the stationery foil.

OCR for page 683

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 11 is a contour plot of the predicted levels of vorticity after the simulation had been running for approximately 900 time steps. We note the presence of both the vorticity injected into the grid by the boundary condition representing the flapping foils and the induced vorticity shed by the stationery foil. The reduced frequency of the unsteady disturbance (in terms of the chord of the stationery foil) is 3.6, corresponding to a disturbance wavelength of somewhat less than twice the chord.
Figure 13. Vector representation of the unsteady velocity field. induced by the flapping foils. The large arrows represent the unit (freeestream) velocity vectors.
Figure 12 presents the resulting prediction of the unsteady lift We note that the amplitude of the unsteady lift force is quite large relative to that which one would predict from classical theory, when only the nominal magnitude of the unsteady vertical velocity component (a purely transverse gust) is considered. We observe, however, that the unsteady velocity field introduced by the flapping foils differs substantially from that of such a transverse gust. This is made apparent in Fig. 13 which shows a vector field describing the unsteady component of velocity at the same time step as the vorticity contours of Fig. 11. The vector field was formed by subtracting the steady solution for the velocity field from the total unsteady solution. Contrast this with Fig. 14 which gives a similar vector field for a transverse gust of k=2.
Figure 14. Vector representation of the unsteady velocity field. of a foil interacting with a transverse gust of k=2.

OCR for page 683

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Vortex Chopping
Figure 15 presents a time history of predicted contours of vorticity as a distributed vortex interacts with a foil. Here the foil is a NACA 0012 thickness form without camber and with angle of attack set to zero. A C-grid with 99×33 grid lines was used for this calculation.
The vortex was generated by perturbing the inflow boundary in a manner similar to that employed for the flapping foil study. We specified an initial vorticity distribution of the form:
(3)
where Г is the total circulation, δ is the core radius, and x0=(x0, y0) is the initial position of the vortex center.
The boundary velocities were perturbed consistent with the potential flow about such a vortex in the absence of the foil, and the nominal position of the vortex was convected, at each time step, with the uniform freestream.
For the results presented in Fig. 15, Г=0.4, δ=0.1 and x0=0. The vortex is seen to intersect with and be split by the foil. The distortion of the vortex is associated with the velocity induced by the foil, which is experiencing lift due to the presence of the velocity field of the vortex.
The results obtained with the present method are seen to agree well with numerical results obtained by Lee and Smith [18], who employed a discrete vortex method to investigate this model problem.
Figure 15. Vorticity contours depicting time history of a distributed vortex interacting with a foil.

OCR for page 683

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
where ej is the jth unit vector. Dennis and Schnabel (19) point out that this is the same as approximating the jth column of the Jacobian by
where a constant value of
was used for this work.
Equation (12) is nothing more than Eq. (30), since the only unknown appearing in this equation is qn+1, that is, in Eq. (31) x=qn+1. Thus, Eq. (31) becomes
N′ (qn+1,m)(qn+1,m+1−qn+1,m)=−N(qm) (33)
where now the aij in Eq. (32) are given by
(34)
It should be noted that the only terms in Eq. (12) that are functions of qn+1 are Δqn and Rn+1. Recall that Rn+1 is defined as
Rn+1=[δi(E−Ev)+δj(F−Fv)]n+1 (35)
or
As mentioned previously, the elements of the Jacobian matrix are calculated by using a finite difference approximation for the derivative. The first two terms of Eq. (13) are the first-order contribution to the higher-order Roe flux vector. It can be seen that the numerical flux at a cell face, say i+1/2 , will be a function of the metrics at i+1/2 and the dependent variables on either side of the cell face. The same statement applies to the viscous flux vector if the cross-derivative terms are not included in its evaluation. Mathematically, for cell face i+1/2, j one has
(36)
and similarly for cell face i,j+1/2
(37)
Now, using the relationships in Eqs. (36) and (37), the left-hand side of Eq. (33) is
(38)
where
and
Note that the first subscript of the Dg (g=e,ev,f or fv) term corresponds to the location of the cell face, and the second subscript corresponds to the location of the dependent variable vector that the numerical flux vector is differentiated with respect to. Also, the right-hand side of Eq. (33) is

OCR for page 683

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
(39)
where
Two items must be addressed pertaining to Eq. (39). First of all, the multiplication of a portion of the equation by the diagonal matrix Ia insures that a divergence free velocity field will be obtained when the iteration in m converges (6). Consequently, the temporal discretization of the continuity equation is taken to be first-order accurate, which leads to the definition of the matrices Ib and Ic in Eqs. (38) and (39). The second item deals with the evaluation of the term , specifically Eq. (11b) must be evaluated as follows to insure that the conservation property at block interfaces is satisfied for blocks in motion (8).
(40)
Solution of Linear Systems
A direct solution of the linear system of equations given by Eqs. (38) and (39) is in general not feasible. One way to circumvent this difficulty is to use a relaxation scheme to solve the linear system at each iteration of Newton's method Using the notation of Reference (18), Eqs. (38) and (39) can be written as
(L+U+B)z=b (41)
where L is a lower block triangular matrix with zeroes on the diagonal, B is a block diagonal matrix, and U is an upper block triangular matrix with zeroes on the diagonal. In Reference (18), the solution of Eq. (41) was computed using the symmetric Gauss-Seidel method A forward pass solution for the change in dependent variable vector z was written as
(L+B)z(1)+Uz(0)=b (42)
where z(0)≡0. Without using z(1) to update b, a backward pass was made using
Lz(1)+(B+U)z(2)=b. (43)
Using Eq. (42), Eq. (43) can be written as
(B+U)z(2)=Bz(1) (44)
or
BΔz =−Uz(2)
where
Δz=z(2)−z(1).
In this work, this procedure has been slightly modified. The forward pass is given by
(L+B)z(1)+Uz(0)=b (45)
where now z(0) ≠0 but is the solution from the previous time step. Once again, z(1) is not used to update b and the backward pass is given by
Lz(1)+(B+U)z(2)=b
or Bz(2)=b−Lz(1)−Uz(2). (46)
All of the solutions contained herein were obtained with three of these symmetric Gauss-Seidel sweeps per Newton iteration.
Turbulence Modeling
The algebraic turbulence model of Baldwin-Lomax (20), which was used in Reference (3) for three-dimensional incompressible flow has been modified for two-dimensional flow. The turbulent wake model used is that of Renze, Buning and Rajagopalan (21). Multiple peaks in the Fmax function are avoided by restricting the search to a specified distance away from vorticity generating surfaces. Also, the method introduced by Chen (22) to calculate wall shear stress and normal distances for nonorthogonal three-dimensional grids has been incorporated into a two-dimensional framework.
Multiblock Approach
An important aspect of computing the flowfield about complex geometries having disparate components with different topologies is grid generation. In this work a domain decomposition technique is used. Once the computation region is decomposed into distinct subregions, each of which

OCR for page 683

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
is gridded independently, the grid topology for each subregion is locally consistent with each component of the geometry. The mechanism by which the subregions are connected leads to a variety of approaches, such as the Chimera grid technique (23), the patch grid technique (24), and block-structured grid technique (25). All three approaches have their relative strengths and weaknesses. A common difficulty to Chimera and patched grid technique is constructing a proper scheme for transferring information between interface boundaries. Nevertheless, both Chimera and patched grid have been successful in modeling flows around multiple bodies in relative motion. One approach to overcome this difficulty (i.e. transfer of information between regions) is the use of block-structured grid techniques. The blocked grid idea is not new. It has been applied to three-dimensional problems by Weatherill and Forsey (26), Belk (27), and Arabshahi (28). The basic idea behind blocked grids is that the whole flowfield between the surfaces of the Configuration and some outer farfield boundary consists of a set of blocks. The union of these blocks fills the entire flowfield without either gaps or overlaps; moreover, the grid has complete continuity (including slope continuity if desired) of grid lines across the block-block boundaries. This feature provides a relatively simple means of communication between neighboring blocks with stationary interfaces in computational space even though the blocks may be in motion in physical space. The approach adopted here with regard to block-block interface communication is a direct extraction-injection procedure which was taken by Belk (27). This means the information (such as q,Δq) from within the domain of one block can be extracted and then injected as phantom data in an adjacent block, thereby eliminating any error due to approximations at block boundaries.
STATIONARY AND DYNAMIC GRIDS
The two major requirements for the grids used in this study are: 1) to simulate the entire experimental configuration, and 2) to mimic the dynamic motion of the flapping foils. The stationary grid (i.e., the flapping foils are fixed) is generated using the algebraic capabilities of EagleView (29). The entire physical domain is divided into four blocks. Each block consists of an H-type grid comprised of 326×65 grid points with a normalized (by the foil chord length) spacing of 1.0×10−5 off solid surfaces. Vertical grid lines are aligned with experimental measurement locations upstream and downstream of the stationary foil. Figure 3 is a
Figure 3. Experimental Locations and Dynamic Grid Boundaries
schematic representation of the experimental measurement locations near the stationary foil.
The efficiency of reconstructing a new grid at each time step to track the movement of the flapping foils is the primary concern for a dynamic grid. To minimize the number of grid points that had to be redistributed, each dynamic portion of the grid is restricted to approximately 101×33 grid points near the corresponding flapping foil. Figure 3 details the extent of this region. The movement of the foils is accomplished by keeping the outer boundary fixed and allowing the foil to oscillate within the confines of these fixed coordinate lines. Weight functions are used to provide a smooth grid distribution between the flapping foils and the outer boundary. These functions, denoted Wi and Wj, are based on an arc length distribution and assumed values from zero to one in each curvilinear direction. The value zero corresponds to no grid motion, while one implies full rotational motion. Thus, the surface of the foil will have a value of one, which decays to zero as the outer boundary is reached. The following combination of weight functions seemed to result in the smoothest grid distribution
(47)
The new position of each point is updated at each time step using the relation
X(new)=(1−Wij) X(fix)+Wij Lcos(α)
Y(new)=(1−Wij) Y(fix)+Wij Lsin(α), (48)
where fix denotes the original position at zero degrees angle of attack, L is the distance from the center of rotation, and α is the angle of attack. Since the grid reconstruction is restricted to a small region and an algebraic method is used to generate the new grid in this region, very little additional computational time is required to update the grid at each time step.
POST-PROCESSING
The computed results had to go through an elaborate post-processing procedure in order to be analyzed in an accurate and meaningful manner. In

OCR for page 683

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
particular the computed results had to be re-referenced to match MIT's reference conditions, and a harmonic analysis procedure had to be used to reduce the large quantity of data obtained from unsteady flows to a manageable set of data that could be easily understood
The measured and computed data are both referenced to ”freestream conditions”. However, evaluation of the measured total pressure coefficient indicates that the reference velocity and pressure , used by MIT, are not measured at the same location in the tunnel. As a result, and are not interchangeable with and , the reference conditions used in the computations. Thus, to compare “apples with apples”, it is necessary to re-reference the computed data to MIT's reference conditions.
The Fourier transform has proven to be an efficient tool for analyzing unsteady data. It works under the premise that any physical process can be described either in the time domain or the frequency domain. The current implementation uses a modified version of a discrete Fourier transform (30) in order to reduce the computing time required to perform the transform. The transform used can be defined as follows
j=1, 2,… N (49)
where
and
k=1, 2,…, kmax.
N=# of time steps per period.
In terms of amplitude and phase, Eq. (49) can be written as
(50)
where
RESULTS
As stated earlier, the computations are being compared with experimental data obtained by MIT (31). Both the steady-state and unsteady computations were carried out using the complete experimental configuration consisting of the water tunnel walls, flapping foils and stationary foil (Fig. 2). The location of the stationary foil coincides with the latest experimental information, specifically, the trailing edge is one centimeter below the centerline of the tunnel and the geometric angle of attack is 1.18 degrees. No slip boundary conditions were used on all solid surfaces, while characteristic variable boundary conditions were used at both the upstream inflow and downstream outflow boundary. Velocities were specified at the inflow while static pressure was set at the outflow. The artificial compressibility factor β was set to 10 for all calculations. All computations were for a Reynolds number of 3.78 million based on the stationary foil chord. The boundary layer on the stationary foil (both the suction and pressure sides) was tripped at a distance of 0.105 chord length from the leading edge, while the boundary layer on the flapping foils and tunnel walls was treated as completely turbulent. A nondimensional grid spacing of 1.0×10−5 off solid surfaces resulted in y+ values of approximately 1.
Steady Results
For the steady-state computations, the flapping foils were held fixed. Since a steady state solution

OCR for page 683

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
was sought, local time stepping was used and the Jacobian matrices were updated every 20 cycles. One Newton iteration was used at each time step as well as three Gauss-Seidel sweeps.
Bounding Box Data
Computed velocities and static pressures were compared with measured data obtained for the bounding box shown in Fig. 3. Data for three ”boxes” was available, but since the boxes were so close together, the data were nearly identical for each box. As a result, the computed results were compared with the data from the middle box. The location of each face relative to the trailing edge of the stationary foil and nondimensionalized with the foil chord is as follows: upstream face, 1.259; downstream face, −0.163; top face, 0.219; and bottom face, −0.153. The static pressures were measured at different locations than the velocities. The upstream static pressures were measured at 0.266 chord length in front of the leading edge of the stationary foil, while the downstream static pressures were measured at 0.156 chord length behind the trailing edge of the stationary foil.
Figure 4 compares the u-component of velocity. The agreement with measured data is good for both the upstream and downstream locations. Figure 5 compares the measured and computed static pressure coefficient at the upstream and downstream faces defined above. The largest discrepancy occurs as one moves from the position of the stationary foil to the upper and lower tunnel walls.
Figure 4. u-Velocity Component
Surface Data
Measured data was available at various locations on the surface of the stationary foil. A comparison of measured and computed surface pressure distributions is shown in Fig. 6. Two representative u-velocity profiles for both the suction and pressure side of the stationary foil are given in Fig. 7. The agreement between experimental and computed profiles on the pressure side, as well as on the suction side upstream of the x/c=0.900 location is good. However, it clearly deteriorates after this station. An argument could be made that this lack of agreement could entirely be attributed to the use of the Baldwin-Lomax turbulence model. However, the inability to match the experimental surface pressure distribution suggests that there may be additional factors. Assuming that the location of the foil in the tunnel and the angle of attack are correct, perhaps there is some slight natural flow angularity in the test section. Another plausible explanation might be that the actual machined geometry differs somewhat from the coordinates given in the experimental reports.
Figure 5. Static Pressure Coefficient
Figure 6. Steady State cp Distribution

OCR for page 683

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 7. Steady State Velocity Profiles
Unsteady Results
The unsteady solution used the steady state solution as an initial condition and the motion of the flapping foils was initiated by boundary conforming dynamic grids that moved in pitch with the flapping foils at a reduced frequency of 3.62 based on semi-chord of the stationary foil. A minimum nondimensional time step of 8.7.×10−4 was used, which corresponds to 1000 time steps per period of motion for the flapping foils. The Jacobian matrix was updated every cycle and three Newton iterations were performed at each time step.
It should be noted that during the course of this work, several different time steps, both larger and smaller, were tried as well as different numbers of Newton iterations. It was found that nearly identical results were obtained for time steps both larger and smaller than the one reported here. However, different numbers of Newton iterations did affect the solution. Newton iterations of one, three and six were used for each time step. Little or no difference in the bounding box results could be seen for any of the different numbers of Newton iterations used. However, a distinct improvement in the viscous (surface) results could be seen between the three Newton iteration solution and the one Newton iteration solution. No significant improvement could be seen between the six and the three Newton iteration solution for either the bounding box or the viscous results.
Although no results were presented here, a grid with nearly twice the resolution in the direction normal to the flappers and foil was also used. The difference in solutions between those obtained from this grid and those presented here was virtually indistinguishable.
Bounding Box Data
The locations for the bounding box used in the steady state analysis are the same for the unsteady analysis, except for the downstream face. The downstream face is located at −0.108 chord length relative to the trailing edge of the stationary foil. For the unsteady case, the upstream and downstream faces of the box were extended into the wake regions. This provided very valuable experimental wake data for comparison with computations.
A comparison of the measured and computed results for the amplitude and phase of the first harmonic of the u-velocity are shown in Fig. 8 for all the faces of the bounding box. The agreement is good, especially on the upstream and downstream faces. It appears that the disturbances generated by the flapping foils are being convected downstream and are arriving at the upstream face with the correct amplitude and phase. Although the agreement at the downstream face does not appear to be quite as good, one can consider the phase at −180 degrees to be +180 degrees and the agreement is quite good.

OCR for page 683

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 8. Unsteady Velocity on Bounding Box-First Harmonic
Surface Data
The mean surface pressure distribution is plotted in Fig. 9, while the amplitude and phase of the unsteady surface pressure distribution is shown in Fig. 10. The discrepancy between the measured and computed mean surface pressure distribution is approximately the same magnitude as the steady state case. The amplitudes of the unsteady surface pressure distribution compare favorably for both the suction and pressure side, while the overall trend of the phase is captured.
Figure 11 contains mean velocity profiles from the suction and pressure side. Once again, the difference between the computed and measured values is about the same as in the steady state case. The amplitude and phase of the first harmonic of the unsteady u-velocity profiles appear in Fig. 12. The computed amplitudes compare favorably with the measured data, while the phase for the pressure side agrees better than that of the suction side. Overall, the comparison of the mean unsteady results showed the same type of discrepancies as the steady results. Figure 13 shows a snapshot of the unsteady velocity distribution.
Figure 9. Mean cp Distribution
Figure 10. Unsteady cp Distribution— First Harmonic

OCR for page 683

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 11. Mean Velocity Profiles

OCR for page 683

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 12. Unsteady Velocity Profiles-First Harmonic
CONCLUSIONS
A two-dimensional multiblock unsteady incompressible Navier-Stokes algorithm based on artificial compressibility has been presented. The unsteady solution was obtained using a dynamic grid that was regenerated at each time step in order to mimic the motion of the flapping foils. For both the steady and unsteady cases, the computed results matched fairly well with the measured bounding box data. However, the computed results for the surface data (in particular, the suction side data) did not match the measured data to the same extent.
Since the solution appears to be grid independent, it seems as though at least two possibilities exist Either the numerical simulation is not being conducted at precisely the same conditions as the experiment or, this is the best that can be expected from the technology in the current algorithm. As to the former, assuming that the angle of attack and location of the foil in the tunnel are accurate, possibilities could include; flow angularity in the test section, alterations in the geometry of the manufactured foil, and the type of flow over the flappers. From a numerical standpoint, the differences could be attributed to the turbulence model used.
Currently, work is being done to adapt the grid to the unsteady wake of the flapping foils in order to minimize the number of grid points required. Future work involves examining alternate turbulence models to see if the discrepancy in the results obtained on the surface can be resolved.
ACKNOWLEDGEMENTS
This work was supported by the Office of Naval Research, grant N00014 –92-J-1060 with James A. Fein as the technical monitor. The authors wish to express their sincere appreciation to the following organizations for providing CRAY-YMP computer resources: U.S. Army Waterway Experiment Station and Cray Research, Incorporated.
Figure 13. Unsteady Velocity Distribution

OCR for page 683

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
BIBLIOGRAPHIC REFERENCES
1. Rice, J.Q., “Investigation of A Two Dimensional Hydrofoil in Steady and Unsteady Flows,” M.S.Thesis, Massachusetts Institute of Technology, June 1991.
2. Delpero, P.M., ”Investigation of Flows around a Two Dimensional Hydrofoil Subject to A High Reduced Frequency Gust Loading,” M.S.Thesis, Massachusetts Institute of Technology, February 1992.
3. Taylor, L.K. and Whitfield, D.L., ”Unsteady Three-Dimensional Incompressible Euler and Navier-Stokes Solver for Stationary and Dynamic Grids,” AIAA-91–1650, June 1991.
4. Roe, P.L., “Approximate Riemann Solvers, Parameter Vectors, and Difference Schemes, ” Journal of Computational Physics. Vol. 43, pp. 357–372, May 1981.
5. Ortega, J.M. and Rheinboldt, W.C., Iterative Solution of Nonlinear Equations in Several Variables. Academic Press, New York, 1970.
6. Pan, D. and Chakravarthy S., “Unified Formulation for Incompressible Flows,” AIAA-89–0122, January 1989.
7. Beam, R.M., and Warming, R.F., “An Implicit Factored Scheme for the Compressible Navier-Stokes Equations, ” AIAA Journal, Vol. 16, No. 4, pp. 393–402, April 1978.
8. Janus, J.M., “Advanced 3-D CFD Algorithm for Turbomachinery,” PhD Dissertation, Mississippi State University, May 1989.
9. Thomas, P.D., and Lombard, C.K., “Geometric Conservation Law and Its Application to How Computations on Moving Grids,” AIAA Journal. Vol. 17, No. 10, 1979.
10. Taylor, L.K., “Unsteady Three-Dimensional Incompressible Algorithm Based on Artificial Compressibility,” PhD Dissertation, Mississippi State University, May 1991.
11. Roe, P.L. and Pike, J., ”Efficient Construction and Utilization of Approximate Riemann Solutions, ” in Computing Methods in Applied Sciences and Engineering, ed. R.Glowinski and J.L.Lions, 6:499–518, Amsterdam: North Holland, 1984.
12. Hsu, C.H., Hartwich, P.M. and Liu, C.H., “Incompressible Navier-Stokes Solutions for a Sharp-Edged Double-Delta Wing,” AIAA -87–0206, January 1987.
13. Rogers, S.E. and Kwak, D., ”Upwind Differencing for the Time-Accurate Incompressible Navier-Stokes Equations,” AIAA Journal, Vol. 28, No. 2, PP.253–262, 1990.
14. Whitfield, D.L., Janus, J.M. and Simpson, L.B., “Implicit Finite Volume High Resolution Wave-Split Scheme for Solving the Unsteady Three-Dimensional Euler and Navier-Stokes Equations on Stationary or Dynamic Grids,” Engineering and Industrial Research Station Report, MSSU-EIRS-ASE-88 –3, Mississippi State University, Mississippi State, February 1988.
15. Sokolnikoff, I.S. and Redheffer, R.M., Mathematics of Physics and Modern Engineering, McGraw-Hill, Inc., New York, NY, 1966.
16. Gatlin, B., “An Implicit, Upwind Method for Obtaining Symbiotic Solutions to the Thin-Layer Navier-Stokes Equations,” PhD Dissertation, Mississippi State University, August 1987.
17. Vanden, K.J., “Direct and Iterative Algorithms for the Three-Dimensional Euler Equations, ” PhD Dissertation, Mississippi State University. December 1992.
18. Whitfield, D.L., and Taylor, L.K., “Discretized Newton-Relaxation Solution of High Resolution Flux-Difference Split Schemes,” AIAA-91–1539, June 1991.
19. Dennis, J.E., Jr. and Schnabel, R.B., Numerical Methods for Unconstrained Optimization and Nonlinear Equations , Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1983.
20. Baldwin, B.S. and Lomax, H., “Thin-Layer Approximation and Algebraic Model for Separated Turbulent Flows,” AIAA-78–257, January 1978.
21. Renze, K.J., Buning, P.G. and Rajagopalan, R.G.,” A Comparative Study of Turbulence Models for Overset Grids,” AIAA-92–0437, January 1992.
22. Chen, J.P., “Unsteady Three-Dimensional Thin-Layer Navier-Stokes Solutions for Turbomachinery In Transonic Flow,” PhD Dissertation, Mississippi State University, December 1991.
23. Benek, J.A., Steger, J.L. and Dougherty, F.C., “A Flexible Grid Embedding Technique with Application to the Euler Equations,” Proceedings of the 6th AIAA Computational Fluid Dynamics Conference, Danvers, MA, 1983.
24. Rai, M.M., “A Conservative Treatment of Zonal Boundaries for Euler Equation Calculation, ” AIAA-84–0164, AIAA 22nd Aerospaces Sciences Meeting, Reno, NV, January 1984.
25. Thompson, J.F., “A Composite Grid Generation Code for General Three-Dimensional Regions, ” AIAA-87–0275, AIAA 25th Aerospace

OCR for page 683

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Sciences Meeting, Reno, NV, January 1987. AIAA Journal, Vol. 26, No. 3, p. 271, March 1988.
26. Weatherill, N.P. and Forsey, C.R., “Grid Generation and Flow Calculations for Complex Aircraft Geometries Using a Multiblock Scheme,” AIAA-84–1665, June 1984.
27. Belk, D.M., “Three-Dimensional Euler Equations Solutions on Dynamic Blocked Grids, ” PhD Dissertation, Mississippi State University, August 1986.
28. Arabshahi, A., “A Dynamic Multiblock Approach to Solving the Unsteady Euler Equations About Complex Configurations,” PhD Dissertation, Mississippi State University, May 1989.
29. Remotique, M.G., Hart, E.T. and Stokes, M.L., ”EAGLEView: A Surface and Grid Generation Program and Its Data Management, ” Software Systems for Surface Modeling and Grid Generation, ed. Robert E.Smith, NASA Conference Publiation 3143, Hampton, VA, April 1992.
30. Singleton, R.C., ”An Algorithm for Computing the Mixed Radix Fast Fourier Transform, ” IEEE Transactions on Audio and Electroacoustics. Vol. AU-17, pp. 93–103, June 1969.
31. Kerwin, J., Keenan, D., Mazel, C., Horwich, and E., Knapp, M., ”MIT/ONR Flapping Foil Experiment, Unsteady Phase,” unpublished data.
DISCUSSION
by Dr. V.C.Patel, The University of Iowa Institute of Hydraulic Research
This question is addressed to this as well as the previous paper. I am very surprised by the rather poor agreement in the velocity profiles with the Baldwin-Lomax model, even in steady flow, especially because this model was optimized for such airfoil flows! Why is that so?
Authors' Reply
Discussions of the tunnel domain with and without Granville's pressure-gradient modifications of the Baldwin-Lomax turbulence model have been included in the paper. It is important to emphasize that Granville's change in the zero-gradient values of Ckleb and Ccp had larger influence in improving the solution than the pressure-gradient corrections to these coefficients.