Direct Numerical Simulation of Unsteady Propulsor Blade Section Forces on a Parallel Processor
W.B.Coney,^{1} S.R.Breit,^{2} and J.R.Webb^{1}
(^{1}BBN Systems and Technologies, USA, ^{2}Kendall 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 twodimensional bladegust interaction problems.
The parallel implementation is discussed in some detail. Timing results obtained on a parallel computer are presented and nearly linear parallel speedups 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 timeharmonic (bladerate and its harmonics) and stochastic unsteady forces produced by interaction with nonuniform inflow are the dominant components of the lowfrequency acoustic signature of a noncavitating 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 bladerate and turbulence ingestion mechanisms which are presently treated separately.
The prediction of the response of a twodimensional blade section to an unsteady gust can be regarded as a preliminary step toward the fully threedimensional 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 bladegust 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 twodimensional bladegust 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.
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 10^{6}),

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 timeaccurate pseudocompressibility (TAP) method with upwindbiased spatial differencing, as described by Rogers & Kwak [7] for our simulations. In the limit where the pseudotime 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 twodimensional model problem, with focus upon one aspect of the bladevortex interaction, is described in Breit, et. al. [9].
The time integration scheme is secondorder accurate and fully implicit, including the boundary conditions. Viscous terms are, for the present, neglected. A sparse, blocktridiagonal, linear system must be solved at each pseudotime 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 completeexchange bandwidth and high randomaccess 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 16byte 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 (i1,j,k), (i,j1,k) and (i,j,k1) 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 speedup of this algorithm. The maximum speedup 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)
For example, on an N×N×N grid, there are 3N2 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 3N2, but the operations can be done in parallel on points which lie on the same plane as defined above.
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/4N^{2}, 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 speedup 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 64^{3}.
The solid line in Fig. 2 represents the speedup for a computation of ideal parallel efficiency, η,
.
(2)
Note in Fig. 2 that nearly ideal speedups are achieved for as many as P=32 processors, with quite good speedup (η=91%) obtained for the maximum number of processors, P=55, utilized in this study.
The falloff 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 dropoff 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 shows some parallel timing results obtained with the present code. Results are for various unsteady simulations performed on the 2D 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 timeshared, multiuser 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 speedups obtained by performing the computations in parallel. The speedup 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 Cgrid with 281×71 grid lines.
Figure 5 gives the steadystate 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, C_{L}=L/(1/2ρV^{2}c), obtained with the two methods are nearly identical, C_{L}=0.570 from the present method, vs. C_{L}=0.576 from the boundary element method.
.
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: onehalf 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 Cgrid 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 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 quasisteady 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
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.
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 nonlinear 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 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 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.
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 Cgrid 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 x_{0}=(x_{0}, y_{0}) 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 x_{0}=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.
CONCLUSIONS
A numerical method for the direct simulation of inviscid unsteady blade section forces has been implemented on a parallel processor. Timing results on the parallel computer showed nearly linear speedups with the number of processors employed. The method was applied to three model problems and the results compared favorably with available solutions.
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.
ACKNOWLEDGMENTS
This work has been supported by the ONR Accelerated Research Initiative on Hydroacoustics of Unsteady Flows through Contract No. N00014– 89C0254. Parallel computer resources utilized in this work were made available by the Lawrence Livermore National Laboratory (LLNL) Massively Parallel Computing Initiative (MPCI). We thank Dr. ChingYeh Hsin and Professor Justin E.Kerwin of the Massachusetts Institute of Technology for providing the boundary element code described in [6].
REFERENCES
1. W.Sears, “Some Aspects of NonStationary Airfoil Theory and its Practical Application, ” Journal of the Aeronautical Sciences. 8(3) 104–108, 1941.
2. H.Naumann & H.Yeh, “Lift and Pressure Fluctuations of a Cambered Airfoil Under Periodic Gusts and Applications in Turbomachinery,” Journal of Engineering for Power, 1–10, January, 1973.
3. H.Atassi, “The Sears Problem for a Lifting Airfoil Revisted—New Results,” Journal of Fluid Mechanics, 141, 109–122, 1984.
4. H.Atassi & T.Akai, “Aerodynamic and Aeroelastic Characteristics of Oscillating Loaded Cascades at Low Mach Number,” Journal of Engineering for Power. 102, 344–351, 1980.
5. H.W.Chiang & S.Fleeter, “Cascade Aeroacoustics Including Steady Loading Effects,” Noise Control and Engineering Journal, 34(2), 61– 72, 1990.
6. C.Y.Hsin, “Development and Analysis of Panel Methods for Propellers in Unsteady Flow,” Ph.D. Thesis, MIT Department of Ocean Engineering, 1990.
7. S.Rogers & D.Kwak, “Upwind Differencing Scheme for the TimeAccurate Incompressible NavierStokes Equations,” AIAA Journal, 28, 253– 262, 1990.
8. M.Athavale & C.Merkle, “An Upwind Differencing Scheme for TimeAccurate Solutions of Unsteady Incompressible Flow,” AIAA Paper 88–3650, 1988.
9. S.Breit, W.Coney, A.Dickinson & J.Webb, “Computing Boundary Forces Due to Unsteady, Inviscid Incompressible Flow,” AIAA Journal, 30(3), 592–600, 1992.
10. S.Breit, W.Coney, J.Webb, “Vortex Chopping Noise in Turbomachinery: Numerical Simulation on a Parallel Processor,” presented at the ONR/DTRC Workshop on Quiet Turbomachinery Design , Annapolis, MD, November, 1991.
11. S.Breit, W.Celmaster, W.Coney, et. al. “The role of Architectural Balance in the Implementation of the NAS Parallel Benchmarks on the BBN TC2000 Computer,” FEDVol. 156, CFD Algorithms and Applications for Parallel Processors, ASME, 1993.
12. S.Breit, “Implementing the Incomplete LU Factorization Method on the BBN TC2000 Parallel Computer,” Proceedings of the Conference on “Parallel Computational Fluid Dynamics.” K.G. Reinsch, et al., ed., Stuttgart, Germany, June 1991.
13. D.Bailey, J.Barton, T.Lasinski and H. Simon, The NAS Parallel Benchmarks. NAS White Paper, 1991.
14. R.Coleman, INMESH: An Interactive Program for Numerical Grid Generation. David W. Taylor Naval Ship Research and Development Center Report DTNSRDC85/054, August 1985.
15. J.Rice, Investigation of a TwoDimensional Hydrofoil in Steady and Unsteady Flows, SM thesis, Massachusetts Institute of Technology, Department of Ocean Engineering , June 1991.
16. P.Delpero, Investigation of Flows around a Two Dimensional Hydrofoil Subject to a High Reduced Frequency Gust Loading, SM thesis, Massachusetts Institute of Technology, Department of Ocean Engineering , February 1992.
17. “MIT/ONR Flapping Foil Workshop”, Carderock Divison/NSWC, Bethesda, MD, March 29–30, 1993.
18. D.Lee and C.Smith, “Distortion of the Vortex Core During Blade/Vortex Interactions,” AIAA 19th Fluid Dynamics, Plasma Physics and Lasers Conference, AIAA87 –1243, 1987.
APPENDIX:
Flapping Foil Boundary Condition
A far field boundary condition model is constructed with potential flow singularities to model the primary elements of the flapping foil experiment. Parameters in the model are estimated using a nonlinear, least squared error approximation to unsteady data and a direct calculation of the steady circulation. All elements of the model are assumed to be two dimensional. The principal elements of the model are:

A single, steady vortex representing the steady lift component on the fixed foil. The strength of the vortex is computed from the measured steady velocity field.

Three unsteady vortices, each representing one of the two moving foils or the unsteady lift component on the fixed foil. Each vortex has a sinusoidally varying circulation distribution in time as in Eq. A1. The magnitude, frequency and phase of the time variation are estimated from the experimental data.

A straight, semiinfinite vortex wake attached to each of the unsteady point vortices. Vorticity is assumed to convect along the length of the wake at a constant rate which is estimated from the experimental data. The local strength of the wake is coupled to the time history of the associated point vortex and the convection velocity and consequantly is sinusoidal in both time and space. Eq. A3 defines the dependence between the point vortex strength and the wake strength. The induced velocity is given by an integral over the length of the wake as shown in Eq. A4.

In order to ensure the velocity is defined at all locations where a computational boundary is likely to be positioned, the unsteady wakes are assumed to have finite, nonzero thickness. If a boundary point falls inside this region, the velocity due to the wake is computed to be the average of the velocity induced at the top and bottom of the region defined by the wake thickness.

Effects due to free vorticity outside of the unsteady wakes and thickness effects of the various airfoils are not included.
The software implementing the boundary condition is provided in two parts. One is a subroutine called by the flow solver. The other is a program used to generate estimates of the parameters in the model from experimental data.
Induced Velocity Formulation
The basic element of the boundary condition model is a point vortex with time varying circulation and an attached wake. The circulation of the point vortex is assumed to be of the form:
Г_{v}(t)=Г_{0}sin(ωt+α) (A1)
and the associated wake will be:
Г_{w}(x,t)=Г′ sin(kx+β) (A2)
where we expect β to be time dependent. The strength of the wake will be related to the strength of the point vortex in the following way. Assume a constant convection velocity for the vorticity in the wake, U. For an infinitesimal change in the strength of the point vortex, dГ_{v}, there is also, at some position Udt away, a filament of strength dГ_{ν}. Since U is constant, we can say:
(A3)
We will be interested in the induced velocity at a particular point in time so we can treat Гw as a function of the distance along the wake. From these assumptions we can see:
From Eq. A2 we can see:
Г′=Г_{0}ω
and
The induced velocity from a point vortex is computed from:
for the coordinate system defined Fig. A1. The vortex is assumed to be located at the origin of the coordinate system and the velocity is computed at point P.
For the wake, the situation is considerably more complicated. For a circulation distribution defined by Eq. A2 and the coordinate system of Fig. A2, the induced velocity is given by:
where and are the unit vectors in the x and y coordinate directions. (The above epression is obtained by integrating the velocity due to a planar vortex sheet once in the direction parallel to the vortex filament we are representing with a point vortex.)
In cases where kx>0 and ky≈0, it is assumed the computation point is inside the viscous region of the wake. For these cases, the velocity is computed at the locations where y=±1/2 (wake thickness) and averaged. The wake thickness is a specified parameter.
Computing the Wake Induced Velocity
Eq. A4 is evaluated numerically over a finite extent of the wake. In order for this approach to be valid, the influence of a semiinfinite wake must approach zero as the point of interest moves far upstream of the start of the wake. With this assumption, we can rewrite Eq. A4 as:
(A5)
where we have assumed kx»ky. We can rewrite this as:
(A6)
and since we also require x<0, the singularity associ ated with is not a difficulty.
We can say:
(A7)
(A8)
for a>0. The sine and cosine integral functions have been used here where:
We are interested in the case where x approaches negative infinity. The integral functions evaluate to:
(A9)
(A10)
where we will take the real part of Eq. A10. Substituting these into Eq. A6 results in:
(A11)
This argument suggests using a finite portion of the wake will eventually converge as longer extents are integrated. No rigorous assesment of the convergence properties of the sequence generated by this process has been performed. Simple numerical experiments suggest it may be necessary to take 1000 or more periods of the wake to obtain accuracy in the range of 1 part in 10^{4}.
Computation of Unsteady Viscous Flow with Application to the MIT Flapping Foil Experiment
E.Paterson and F.Stern
(University of Iowa, USA)
ABSTRACT
A timeaccurate unsteady viscousflow method is validated through calculations and comparisons with the Massachusetts Institute of Technology flappingfoil experiment. Solutions are obtained using a small domain surrounding the foil, a tunnel domain that included the foil and the tunnel walls, and a complete domain that included the foil and both the tunnel walls and the upstream flappers. In the latter case, the CHIMERA overlaidgrid method was used. The solutions give similar overall agreement with the data for both steady and unsteady flow, which demonstrates that such problems can be handled with a variety of formulations, although the boundary data, cpu time, and storage requirements are different. The physics are complex with analogy to Stokes layers and explicated through analysis of the axial pressure gradient, which exhibits upstream and downstream traveling waves over the foil and in the near wake and in the intermediate wake, respectively, due to nonlinearities induced by the convective acceleration and steady/unsteady interactions. The nature of the unsteady displacement thickness suggests viscousinviscid interaction as a possible mechanism for the axial pressuregradient response.
NOMENCLATURE
c 
=wave speed (=λ/T) 
C_{p} 
=pressure coefficient 
L 
=characteristic (foil) length 
p 
=pressure 
Re 
=Reynolds number (=Uo_{L}/v) 
t 
=time 
U_{i} 
=ensemble averaged velocity components in Cartesian coordinates (= U,V,W) 
=Reynolds stress tensor 

x^{i} 
=Cartesian coordinates (=x,y,z) 
y^{+} 
=wall coordinate (=Uτy/v) 
δ 
=boundarylayer thickness 
γ 
=phase angle 
ν 
=kinematic viscosity 
ν_{t} 
=eddy viscosity 
ξ 
=frequency parameter (=ωL/U_{0}) 
ξ^{i} 
=curvilinear coordinates (=ξ,η,ζ) 
τ_{w} 
=wallshear stress 
=transport quantities (U,V,W) 

ω 
=circular frequency (=2πf) 
INTRODUCTION
Although most practical engineering flows are inherently unsteady and viscous, the physical understanding is in its infancy in comparison to steady flow, due to the complexity of performing experiments and computations even for simplified configurations, as evidenced by the sparse number of studies documented in the literature for unsteady vs. steady flow (e.g., the proceedings of this conference). Unsteady flows can be conveniently grouped into three categories: (1) natural, (2) forced, and (3) interacting. For all three, the predominant number of studies are experimental vs. computational, as indicated by the example references given.
Naturally unsteady flow refers to conditions under which the unsteadiness arises from flow separation and associated vortex wakes or jet and shearlayer instabilities. Most information is for twodimensional flows, which includes the Reynolds number (Re) dependence of the separation points, vortexshedding frequencies [i.e., Strouhal number(St)], and wakepatterns for a variety of bluff bodies both with sharp (e.g., flat plates at incidence) and smooth (e.g., circular cylinders) edges [1,2,3]. The situation is similar for threedimensional flows, but the number of studies is drastically reduced [4,5].
Forced unsteady flow refers to conditions under which the unsteadiness arises from body motions or moving boundaries (e.g., rotortype flows and control surfaces) or oscillating streams (e.g., boundary layers with oscillating external flows) such that there is an identifiable predominant forcing frequency. With respect to the former, relatively few studies have been conducted for realistic geometries, due to the complexity of such configurations and the geometric dependence of the resulting data. However, a number of studies have examined the idealized geometry of a pitching foil where the interest has ranged from inviscid flutter analysis [6] to dynamic stall [7–9]. The situation is similar for oscillating external flows where model problems such as flatplate boundary layers (i.e., Stokeslayer overshoots, phase angles, and streaming) both for laminar and turbulent flow [10] and foils embedded in transverse gusts [11] have been studied.
Interacting unsteady flow refers to conditions under which both natural and forced unsteadiness are present and interact. Obviously, in this case, the least information is available, although this may, in fact, be the most important category with regard to engineering design improvement. The few studies that have been conducted were concerned with wake and separation control [12,13] and blade/vortex [14,15] and rotor/stator [16] interactions. In the former cases, it has been demonstrated that wakes with harmonic disturbances behave like nonlinear oscillators and are characterized by synchronization (i.e., lockin) vs. non lockin conditions depending on the ratio of the frequency of the disturbance f_{e} to St (i.e., stateselection diagram or resonant horns of amplitude vs. f_{e}/St), which has been further explicated through stability analysis [ 17], and similarly the ability to control separation through disturbances at certain f_{e}. In the latter cases, the studies are somewhat preliminary in that, thus far, the blade and stator flows, respectively, are usually steady.
The present work is part of a larger project at the Iowa Institute of Hydraulic Research (IIHR) specifically concerning the development of computational fluid dynamics (cfd) methods with application to marine propulsors. Marine propulsors are unique in comparison to related applications (e.g., turbomachinery, turboprops, and rotors) in that they operate in the thick hull boundary layer and/or appendage wakes such that complex naturally (hull flow) and forced (propulsor flow) unsteady interactions occur. Initially, cfd methods were applied to propeller/hull interaction utilizing an interactive approach, i.e., a bodyforce propeller representation in the viscousflow method obtained interactively using a vortexlattice propellerperformance method for specified inflow [18]. Subsequently, cfd methods were applied directly to calculating the propellerblade and wake flows both for idealized [19] and practical turboprop and marine propulsors [20]. Currently, efforts are being directed towards extending the latter methods for unsteady flow with consideration to issues of efficient largescale timeaccurate computations for fixed and movingboundary problems [21, 22].
The present paper describes the initial efforts in extensions for unsteadyflow calculations and validation through calculations and comparisons with the Massachusetts Institute of Technology (MIT) Marine Hydrodynamics Laboratory flappingfoil experiment (ffx). ffx represents a twodimensional simulation of propellerblade and wake flow, i.e., a blade section (foil) embedded in vertical and horizontal gusts generated by upstream pitching foils (flappers) (figure 1). ffx was conducted such that calculations could be performed either for a small domain with given boundary data (sd), for the tunnel domain with a specified inflow (td), or for the complete domain, including the upstream flappers (cd). Associated with ffx was the 29–30 March 1993 Office of Naval Research (ONR)/MIT UnsteadyFlow Workshop in which various groups of researchers submitted blind computations (i.e., with only the experimental conditions and boundary data) for either sd, td, or cd and met at the workshop for the comparisons with the data and discussion. At the workshop, IIHR submitted blind computations for sd. Herein, these results will be presented along with more recent results for sd, td, and cd. In the latter case, the CHIMERA overlaidgrid method was used, which has also been under investigation for other fixed and movingboundary problems [22]. The recent results were presented at the 13 September 1993 post workshop meeting.
In the following, the computational method and ffx are described and the computational conditions, grids, and uncertainty are given. Then, the results are discussed for both steady and unsteady flow, including comparisons with the data, discussion of the ONR/MIT workshop and post workshop meetings, and analysis and analogy to Stokes layers. Lastly, some concluding remarks are made, including future work.
COMPUTATIONAL METHOD
The computational method is based on extensions of [19] for timeaccurate unsteady
flow calculations. Although the basic viscousflow method [23] utilized in [19,20] was developed with the capability of unsteadyflow calculations, only very limited such calculations had previously been performed. In the following, an overview is given of the equations and coordinate system, turbulence model, discretization and velocitypressure coupling, solution domains and boundary conditions, and grid generation. The modifications required for timeaccurate computations are also discussed. The method is applicable to unsteady threedimensional flow; however, the current presentation is for unsteady twodimensional flow in view of ffx. Further details concerning the basic viscousflow method are provided in [23] and associated references.
Equations and Coordinate System
The continuity and Reynoldsaveraged Navier Stokes (RaNS) equations for unsteady incompressible flow are written in the physical domain using Cartesian coordinates and tensor notation with x positive downstream, y normal to x, and the origin at the foil leading edge (figure 1). In nondimensional form the equations are
(1)
(2)
with and U_{i}=(U,V,W) and u_{i}= (u,v,w) are the components of the ensembleaveraged velocity and turbulent fluctuations, respectively, and p is the pressure. All variables are nondimensionalized using U_{o}, L, and ρ.
The Reynolds stress is related to the mean rate of strain by an isotropic eddy viscosity v_{t}
(3)
upon which (2) becomes
(4)
where . (5)
Transformation of (1–5) into nonorthogonal curvilinear coordinates is such that the computational domain forms a simple rectangular parallelepiped with equal grid spacing. The transformation is a partial one since it involves the coordinates only and not the velocity components. The transformation is accomplished through the use of the expression for the divergence and “chainrule” definitions of the gradient and Laplacian operators, which relate the orthogonal curvilinear coordinates x^{i}=(x,y,z) to the nonorthogonal curvilinear coordinates ξ^{ι}= (ξ,η,ζ). In this manner, (1–5) can be rewritten as follows
(6)
(7)
with
(8)
and where is nonzero only for movinggrid applications.
Turbulence Model
Closure is attained through a quasisteady application of the BaldwinLomax turbulence model with modifications to account both for the effects of wake asymmetry [24,25] and axial pressure gradients [26]. The former was used to obtain the solutions for all three domains (sd, td, and cd), whereas the latter was only used for td (denoted td_{dp/dx}).
For the standard BaldwinLomax model, in the nearwall layer
(9)
where the length scale is based on the mixing length
(10)
and the vorticity ω is used to determine the velocity scale. In the outerlayer
(11)
where
(12)
(13)
(14)
(15)
F_{Kleb} is the Klebanoff intermittency function. In F_{wk,} the first and second expressions are effective in attached and separated boundary layers (and wakes), respectively. The velocity scale F_{max}, is determined from (14) and y_{max} is the corresponding ycoordinate. The empirical constants are: κ=0.4, k=0.0168, C_{cp}=1.6, C_{Kleb}=0.3, C_{wk}=0.25, and A^{+}=26.
The standard BaldwinLomax model is not suitable for asymmetric wakes where two different values for y_{max} and F_{max} result from the merging of the suction and pressureside boundary layers. Therefore, the eddy viscosity on the farwake centerline is rendered singlevalued by taking the maximum value from both sides of the wake.
v_{t,fwk}=kC_{cp}F_{Kleb} max(F_{wk,top},F_{wk,bottom}) (16)
The eddyviscosity distribution is asymmetric due to the evaluation of F_{kleb} on each side of the wake. For both F_{Kleb} and F_{wk,}y is measured from the point of maximum wake deficit. Finally, since (11) is only valid in the far wake, a blending function is used to provide a smooth transition from the boundary layer to the far wake
(17)
where v_{t,te} is the distribution at the trailing edge of the foil.
The effects of axial pressure gradient were investigated through the use of Granville's modified coefficients C_{cp} and C_{Kleb} such that they are functions of a Clausertype pressuregradient parameter
Note that the zero gradient values C_{cp}=1.2 and C_{kleb}=0.646 are significantly different from the standard BaldwinLomax values.
Discretization and VelocityPressure Coupling
The finiteanalytic method is used to reduce (7) to algebraic form. Each local rectangular numerical element is linearized by evaluating the coefficients and source functions at the interior node P and transformed again into a normalized form by a simple coordinate stretching. In addition, the temporal derivative is expressed as a secondorder backward finite difference. An analytical solution is obtained by separation of variables with specified boundary functions. As a result, a 9point finiteanalytic formula is obtained
(18)
The resulting computational stencil includes all eight neighboring nodal values and the values at the two previous timesteps. Since (18) is implicit, both in space and time, its assembly for all elements results in a set of simultaneous algebraic equations that is solved using a lineADI approach with underrelaxation.
The coupling of the velocity and pressure fields in an efficient timeaccurate fashion required changing from a SIMPLER algorithm, which requires global iterations for time accuracy, to a splitoperator approach based on the PISO algorithm. The first step is to implicitly solve the momentum equations given the pressure field from the previous time step using (18). Next, the implicit pressure equation and the explicit velocity correction are solved iteratively until continuity is satisfied. In both steps, the finiteanalytic coefficients C_{nb} are updated in each iteration, thus retaining the nonlinear nature of the RaNS equations.
Solution Domains and Boundary Conditions
In terms of the notation of figure 1, the boundary conditions are as follows.
For sd: on the top and bottom surfaces S_{sd1},
(19)
where irrotationality has been assumed so that the velocity derivatives in the evaluation of the pressure can be prescribed by the data and the upstream lower corner of the sd boundary (see figure 1) is fixed at an instantaneous reference pressure defined as the volume average over the computational domain [27]
(20)
on the inlet S_{sd2},
(21)
and on the exit plane S_{sd3},
(22)
where U_{e} and V_{e} are specified from the data.
For td: on the tunnel walls S_{td1},
(23)
on the inlet S_{td2},
(24)
where U_{e} and V_{e} are a combination of the data (figure 3) and an inviscid model of an internal progressive wave (Appendix A), and on the exit plane S_{td3},
(25)
For cd: on the tunnel walls S_{cd1}and on the exit plane S_{cd3}, (23) and (25) are applied, respectively, on the inlet, S_{cd2},
(26)
i.e., uniform flow is assumed, and on the flapper surfaces S_{cdf}.
(27)
where α=α_{1} sin(ξτ) is the flapper angle, (x_{r},y_{r}) is the center of rotation, τ is the nondimensional time, and ξ (=ωL/U_{o}) is the frequency parameter (or St), which is twice the reduced frequency k (=1/2ξ).
Grid Generation
Three gridgeneration techniques were used depending on the domain and grid topology. EAGLE [28] was used to generate cgrids for both the foil in sd (figure 2) and the flappers in cd (figure 4). Algebraic hgrids were used for td and cd (figures 3 and 4) with hyperbolic tangent stretching functions and transfinite interpolation such that the grid spacing was controlled along the foil surface, leading and trailing edges, and the flapper wakes. CHIMERA domain decomposition [29,30] was used for cd to resolve the oscillating flappers by decomposing cd into three separate, but coupled domains. Gridtogrid communication was facilitated through bilinear interpolation of the flow variables (U,V,p) from the foil grid to the flapper outer boundaries and from the flapper grids to the foil hole boundaries. Holes in the foil grid are created by excluding points that lay within specified boundaries on the flapper grids. Typically, these boundaries are specified such that a three or four point overlap region exists. The interpolation coefficients and holes and hole boundaries were determined a priori for each time step. Figure 4 shows global and detailed views of the overlaid grid system.
The use of oscillating, overlaid grids required modifications of the tridiagonal solver
to incorporate the hole boundary conditions and to solve across hole regions
a_{i}=a_{i}• IBLANK_{i}
b_{i}=b_{i}• IBLANK_{i}+(1−IBLANK_{i})
c_{i}=c_{i} • IBLANK_{i}
d_{i}=d_{i} • IBLANK_{i}+(1−IBLANK_{i})•φ_{i} (28)
where IBLANK is 0 for a hole point and 1 for a field point, a_{i}, b_{i}, c_{i}, and d_{i} correspond to the lower, center, upper diagonals, and right hand side of the tridiagonal matrix, respectively, and φ_{i} corresponds to the flow variable value inside or on the hole boundary.
MIT FLAPPINGFOIL EXPERIMENT
ffx represents a twodimensional simulation of propellerblade and wake flow, i.e., a foil embedded in vertical and horizontal gusts. The gust (i.e., travelingwave) models the flow experienced by a propellerblade section as it rotates through the nonuniform hull and/or appendage wake. The purpose of the experiment was to provide detailed measurements for validation of unsteady cfd methods and to provide data for the determination of appropriate Kutta conditions for unsteady potentialflow methods. Details of the experimental objectives, apparatus and procedures are documented in [31–33].
Experimental Geometry, Information, and Uncertainty
The experiments were performed in the closedloop MIT VariablePressure Water Tunnel. The test section is 53 inches long and 20 inches square with a freestream turbulence level of about 1%. The foil is a modified NACA 66 with a chord and span length of 18 and 20 inches, respectively, and was fixed at an angle of attack of 1.18–1.34º. Transition was set at x=0.105 on both sides using 0.050 inch diameter epoxy disks 0.008 inch high separated by 0.050 inch. The flappers are NACA 0025 with a chord and span length of 3 and 19.5 inches, respectively, and are driven by a variablespeed motor, which was designed to operate at reduced frequencies of up to 10. The flappers oscillated with an amplitude of 6º and a frequency of 16 Hz. The corresponding Re based on the foil chord length and reference velocity 20.94 ft/sec is 3.78 ×10^{6} and the reduced frequency is k=3.6
Velocity and surfacepressure measurements were made using a twocomponent laser Doppler velocimeter and miniature pressure transducers, respectively. For the unsteady measurements, the value at each time step is the ensemble, or phase, average over 250 periods and had a temporal discretization of 180 points over the period.
Prior to the workshop, steady and unsteady (U,V) data was provided on sd boundaries (figure 2) and a large portion of td inlet (figure 3) along with steady boundarylayer profiles. After the workshop, unsteady boundarylayer profiles and pressure data was provided at the locations shown in figure 1.
Experimental uncertainty is due to instrument accuracy, geometric uncertainty, and data processing. The velocity data on sd boundaries has a standard deviation of 0.029 in the streamwise direction and 0.016 in the vertical direction. Standard deviations for the boundarylayer profiles vary from about 0.15 close to the foil surface to about 0.025 at the edge of the boundary layer. The unsteady pressure measurements had a standard deviation of 0.02 based upon absolute pressure. There was some uncertainty in the foil angle of attack due to the difficulty of determining the precise location of the foil in the tunnel. The data was normalized by the velocity at the exit of the contraction measured without the flappers and foil, i.e., 20.94 ft/sec. This presents no problems for sd and td comparisons since they use the data for boundary conditions; however, cd was normalized based on a unit reference inlet value such that the presence of the flappers and foil are included, which results in a 3% lower velocity on S _{td2} in comparison to the data. A similar difference was confirmed by MIT through inviscidflow predictions.
COMPUTATIONAL GRIDS, CONDITIONS, AND UNCERTAINTY
The grids for sd, td, and cd are shown in figures 2–4, respectively, sd and td grids have 181×80 (14,480) and 180×179 (32,330) points, respectively, cd has 240×179 (42,960) points and 71×40 (2840) points for the foil and flappers, respectively, for a total of 48,640. For each domain, approximately 40 points span the boundarylayer and the nearwall grid spacing was set such that the first grid point was located at y^{+}≈0.1.
The experimental values of Re and ξ were used in the calculations. The calculations were made at an angle of attack of 1.34º. An additional td calculation was made for angle of attack 1.18º, which indicated a 2% change in the steady lift. The unsteady solution was unchanged. The solution submitted for the workshop (denoted sd_{orig}) used the data, which was smoothed using a cubic spline and
interpolated both in time and space along sd boundaries using a biquadratic polynomial, to evaluate the pressure from (19). However, the coarseness of the data introduced errors in the evaluation of (19) and lead to erroneous results, especially for pressure. Subsequently, a revised sd (denoted sd) used td to provide finer resolution data for the boundary conditions. For td inlet boundary S_{td2}, smoothed and interpolated data provided approximately 60% of the tunnel inflow area and, as previously mentioned, the remaining 40% was specified using a potentialflow approximation. For cd, the boundaryconditions were wellposed and required no special treatment. However, no turbulence model was used in the flapper domains. Lastly, for each domain, transition was fixed by forcing the eddy viscosity upstream of the boundarylayer trip to zero.
Timestep sensitivity studies were conducted for sd which showed that approximately 50 time steps/period provided a timestep independent solution. For sd and td, the time step was 1/50 of a period (Δt=0.01745) and for cd, the time step was 1/64 of the period (Δt=0.01364). Grid sensitivity studies were conducted for steady flow from which the boundarylayer and nearwall grid spacing described earlier was determined. Resolution of the point distribution in the wake of the flappers for both td and cd was not studied.
The steadyflow solution was used as the starting point for each calculation and typically 2 periods for sd and 4 periods for td and cd were required to attain a periodic solution. Cpu time per period on a CRAY YMP was 15 minutes, 45 minutes, and 2 hours for sd, td, and cd, respectively, and storage requirements ranged from 4 MW for sd to 8 MW for cd. Typically, 25 and 50 iterations were required for the implicit momentum and the pressurevelocity correction steps, respectively, to satisfy the convergence criterion of R<0.0001, where the residual was defined as
(29)
and i, it, and imax are the gridpoint and iteration indices and total number of grid points, respectively.
Uncertainty in the solutions is due to the boundary conditions and truncation errors. For sd and td, the process of using data for boundary conditions introduces errors from smoothing and interpolating, experimental uncertainties, and insufficient spatial resolution. The numerical method is secondorder in time and approximately first order in space. Therefore, given Δt≈0.01, the temporal truncation errors are O(10^{−4}).
COMPARISON OF SOLUTIONS AND DATA
sd, sd_{orig}, td, td_{dp/dx}, and cd results are compared with each other and the data. Comments are made concerning the overall results at the ONR/MIT UnsteadyFlow Workshop and post workshop meetings. The velocity, pressure, shearstress, and force harmonic amplitude and phase are Fourier analyzed
(30)
where γ_{n} is with respect to the flapper angle.
Steady Flow (Zeroth Harmonic)
Figure 5 shows td and cd (U_{0},V_{0}) on S_{td2} and S_{sd1} top and bottom, including data. On S_{td2}, cd U_{0} shows smaller magnitude and flapper wakes with larger width and deficit than the data and V_{0} in agreement. On S_{sd1}, td and data U_{0} agree, whereas cd agrees in shape, but under predicts the magnitude, td and cd V_{0} are similar, but show an axial decay in comparison to the data.
The surface pressure is shown in figure 6. sd and td are close to the data on the pressure side, but over predict on the suction side, cd is similar on the pressure side and closer to the data on the suction side. The latter is consistent with the lower inlet velocity. Over the adverse pressuregradient region on the suction side, td_{dp/dx} shows improvement (i.e., closer agreement with the data). The corresponding lift for sd, td, cd, sd_{orig}, and td_{dp/dx} are 0.65, 0.65, 0.60, 0.54, and 0.63.
The wallshear stress is shown in figure 7. The fixed transition along with a small region of separation on the last 5% on the suction side are displayed, sd and td agree, but cd shows lower values.
The steady velocity profiles are shown in figure 8. The solutions agree with the data, except near the trailing edge, where td_{dp/dx} shows improvement, especially on the suction side. The change in the zerogradient values of C_{kleb} and C_{cp} had a larger influence than the corrections, cd shows a consistent 3% under prediction outside the boundary layer Wallcoordinate profiles of the data [34] indicate extensive regions of transitional flow on both the pressure (x< 0.76) and suction (x<0.61) sides, (i.e., unexpected departures for the data from the logarithmicoverlap law) and boundary effects on the nearwall data, which may partly explain the differences between the solutions and data. However, the level of agreement for the solutions without and with the pressuregradient modifications appears to be a general assessment of the current capabilities of isotropic turbulence models since it is consistent with the overall results reported in the literature.
Unsteady Flow
Figure 5 includes td and cd firstharmonic amplitude (U_{1},V_{1}) and phase (γ_{U1}, γ_{V1}). On S_{td2}, cd agrees in amplitude shape and phase, but under predicts the former magnitude. Also shown are td secondharmonic amplitudes (U_{2}, V_{2}), which are large, particularly near the flapper wakes. On S_{sd1} top, td U_{1} agrees in magnitude, but lacks the spatial oscillations in the data, cd U_{1} is under predicted, but displays similar spatial oscillations as the data, td V_{1} agrees in magnitude and shape and cd V_{1} in magnitude with the data. On S_{sd1} bottom, both td and cd U_{1} show agreement in magnitude, but under predict the spatial oscillations, td and cd V_{1} show an axial decay in comparison to the data. The phases are in agreement with the data and show a downstream traveling wave [in the Fourier analysis, the argument of the trigonometric functions (ξt+γ) is equivalent to characteristic lines ξ(t±x/c), which corresponds to downstream(γ=−ξx/c) or upstream(γ= +ξx/c) traveling waves].
Figures 9 and 10 show the lift and drag time histories. The magnitudes are similar; however, sd, td, and td_{dp/dx} indicate a second harmonic, whereas cd is primarily first harmonic, and sd_{orig} shows large oscillations. There are larger differences for ΔC_{D} than ΔC_{l},
Figure 11 shows the surfacepressure firstand secondharmonic amplitude and phase. In comparison to the data, sd, td, and td_{dp/dx} show amplitudes of similar magnitude. The shapes and zero values at certain x/L are different and unconfirmed by the limited data, td_{dp/dx} shows improvement. The large second harmonic is consistent with the lift and drag and was attributed by MIT to an excited mode of tunnel vibration [33]. Recall the large td secondharmonic amplitudes on S_{td2}. cd firstharmonic amplitude is larger than the others and cd secondharmonic amplitude, the latter of which is similar to the others. This trend is consistent with the lift and drag. On both sides, cd and td_{dp/dx} indicate increasing phase (i.e., upstream traveling wave), whereas sd and td indicate nearly constant values with 180º changes at the zerovalue locations. The limited data indicate constant phase with no particular agreement with the solutions.
Figure 12 shows the wallshear stress firstand secondharmonic amplitude and phase. The solutions show amplitudes of similar magnitude for both harmonics. The firstharmonic phases indicate: on the pressure side, decreasing values (i.e., downstream traveling wave) on the forebody and upstream traveling wave on the afterbody for sd, td, and td_{dp/dx} and vice versa for cd; and on the suction side, downstream traveling wave for sd and td, cd similar to the pressure side, and upstream traveling wave for td_{dp/dx}. The secondharmonic phases indicate: on the pressure side, nearly constant values with 180º changes at the zerovalue locations for sd, td, and cd and upstream traveling wave for td_{dp/dx}; and on the suction side, nearly constant values with 180º changes at the zerovalue locations for sd, td, and td_{dp}/_{dx} and upstream taveling wave for cd.
The velocity firstharmonic amplitude and phase (U_{1}, γ_{U1}) profiles are also shown in figure 8. Both the solutions and data indicate, for x≤0.784 (pressure side) and ≤ .9 (suction side), amplitudes with small overshoots and phases with increasing lags with x and smooth transition across the boundary layer, whereas for x greater than these values, amplitudes and phases with a twolayer structure: relatively constant values of different magnitude for the inner and outer flows with a zero amplitude and 180º phase shift at the inner/outer boundary (i.e., boundarylayer thickness). The magnitudes are close to the data and td_{dp/dx} shows improvement, especially near the trailing edge. In the wake, the solutions show large differences, especially near the wake centerline. sd, sd_{orig}, and td have a smaller amplitude near the centerline and a larger amplitude outside the boundary layer in
comparison to cd and and td_{dp/dx}, which show large overshoots. Initially, the phase shows a continuation of the boundarylayer response and then by x=1.1, leads for sd, td, and td_{dp/dx} and lags for cd and sd_{orig} on the pressure side and leads on the suction side.
Figure 13 shows axialvelocity contours for sd, td, and cd at t=0. The difference in the domain size, consistency between the solutions, and wavy nature of the flow are evident, cd shows continuity across the overlaidgrid region and flapper wakes that are wider than td, which indicates almost zero width.
The solutions are similar in showing overall good agreement with the data, except for sd_{orig} unsteady lift and drag; however, detailed differences are evident between the solutions and the data. Since td_{dp/dx} showed improvement, it alone will be used for further discussion.
Figure 14 shows the perturbationvelocity (i.e., difference between steady and unsteady) timehistory contours, which vividly exhibit the features described with regard to the firstharmonic velocity profiles and data. td_{dp/dx} trends are smoother, due to the coarseness of the data (13 points across the boundary layer).
Figure 15 shows the perturbation velocity vectors and particle traces at various time steps. The interaction between the flapperinduced vortices and foil is complex: distortion and increased speed on the suction (c=1.26) vs. the pressure (c=1.05) side such that 5 vs. 4 vortices are observed; and secondary counterrotating vortices near the trailing edge and in the wake. The latter directly correlate with and explicate the region where the firstharmonic velocity profiles displayed the twolayer structure.
At the ONR/MIT UnsteadyFlow Workshop [35] and post workshop meetings, ten groups from naval laboratories (DTMB & NRL) and universities (Iowa, Cincinnati, Minnesota, Mississippi State, Penn State, VPI) submitted solutions using a variety of methods (inviscid, boundarylayer, and RaNS) and domains (sd, td, and cd). There were differences/difficulties between/with the solutions for both steady and unsteady flow: the groups using sd had similar problems; the steady lift ranged from 0.54 to 0.70; inability to predict the details of the steady separation; and the unsteadyflow trends varied considerably. However, three groups, Iowa, MSU [36], and PSU [37] had the best results of a similar encouraging level of agreement with the data.
ANALYSIS AND ANALOGY TO STOKES LAYERS
The results are further analyzed through discussion of the relationship between the flow pattern and axial pressure gradient. The analogy between the present flow and Stokes layers for other idealized geometries is discussed.
Figure 16 shows the perturbation axial pressuregradient contours at the same time steps as figure 15. Comparison of the figures indicates the direct correspondence between the flow pattern and axial pressure gradient, i.e., the flow directions are consistent with the regions of favorable and adverse gradients. The complexity of the wake is evident, including higher harmonics and upstream and downstream traveling waves over the foil and in the near wake and in the intermediate wake, respectively. The upstream traveling waves are consistent with the surfacepressure firstharmonic phase.
The nature of this correspondence can be explicated by solving the axial Euler equation [i.e., (2) without the viscous and Reynoldsstress terms] for the axial pressure gradient and evaluating the relative importance of the various terms using Fourier analysis
(31)
and td_{dp/dx}. Corresponding to orders U_{0,}U_{1,} and U_{2}, 2+11+19=32 terms arise, respectively, which were evaluated along lines j=40 and 185 outside the boundary layer at y^{+} ≈10^{4} (figure 17a ). In the latter cases, three terms are dominant, one from the temporal and two from the convective acceleration, such that
(32)
where
(33)
and is with respect to . Three cases are of particular interest: (1) U_{n}≠ U_{n}(x) (i.e., U_{x,n}=0) such that A=±∞ and = ±π/2, where+corresponds to a lead (i.e., temporal wave with c → ∞) and−to a lag (i.e., spatial wave with c → 0); and (2) U_{n}≈ax (i.e., U_{x},n=a) such that A and ; and (3) U_{n}≈sin(nξxB/c) (i.e., U_{x,}_{n}≈cos(nξxB/c)) such that A≈tan(nξxB/c) and where B=1 (i.e., U_{n} primarily first harmonic) corresponds to a temporal wave and B≈2 (i.e., U_{n} primarily second harmonic) an upstream traveling wave.
Figure 17b shows the velocity first and secondharmonic amplitudes along j=40 and 185. The former is primarily second harmonic on the foil (i.e., case 3) and linearly decreasing in the wake (i.e., case 2), which is consistent with figure 5. Figure 17c indicates B=1.2 and 1.8 on the pressure and suction sides, respectively, i.e., upstream traveling waves, which is consistent with figures 11 and 16. The latter is primarily third and second harmonic on the pressure and suction sides, respectively, on the foil and nearly zero in the wake.
Figure 17d shows the displacementthickness firstharmonic amplitude. Similar spatial oscillations as U_{1} are observed, but with relatively larger values in the near wake. The maximum value coincides with the location where changes from upstream vs. zero phase, which suggests viscousinviscid interaction as the mechanism for the pressuregradient response.
Figure 18 shows the velocity, pressure, and pressuregradient firstharmonic phase for the outer (y^{+}≈10^{4}), overlap (y^{+}≈200), and sublayer (y^{+}≈0.1) regions. The pressuregradient phase was obtained by differentiation of the pressure and, additionally, the outer region also includes the Eulerequation solution (33). The pressure and pressure gradient indicate upstream traveling waves. At the outer region, close agreement is shown between the two pressure gradients, which validates the Eulerequation analysis. On the foil in the outer and overlap regions, the velocity shows a downstream traveling wave, whereas in the sublayer region it follows the pressure gradient, which is consistent with the wallshear stress. In the wake, the velocity shows a downstream traveling wave.
Figure 19 shows the velocity first harmonic amplitude and phase in wallcoordinates. The largest overshoots occur near y^{+} ≈1000, where the phase abruptly changes from the outer to inner values. The amplitudes display a doublepeak across the boundary layer and the phase is constant in the sublayer. The large change in phase shown in the overlap region (10<y^{+}< 1000) correspond to the lags shown in the velocity profiles. The role of the turbulence model vs. physics for the sublayer region is unknown, as no data is available.
In conclusion, ffx physics are complex with analogy to Stokes layers (i.e., overshoots and phase shifts) for temporal, spatial, and traveling horizontalwave external flows [10,38,39]. The latter study has shown the differences/similarities between these flows and that the response for the latter two is extreme due to nonlinearities induced by the convective acceleration such that further study is required, especially to document the nature of the interaction between the mean and turbulent motions. The present work supports this conclusion and, moreover, displays the additional complexity for combined vertical and horizontal waves.
CONCLUDING REMARKS
A timeaccurate unsteady viscousflow method has been validated through calculations and comparisons with the Massachusetts Institute of Technology flappingfoil experiment. Solutions were obtained using a small domain surrounding the foil, a tunnel domain that included the foil and the tunnel walls, and a complete domain that included the foil and both the tunnel walls and the upstream flappers. The latter case involved the use of the CHIMERA overlaidgrid method. The three solutions gave similar overall agreement for both steady and unsteady flow, which demonstrates that such problems can be handled with a variety of formulations, although the boundary data, CPU time, and storage requirements are quite different. The physics are complex with analogy to Stokes layers and are explicated through analysis of the axial pressure gradient, which exhibits upstream and downstream traveling waves over the foil and in the near wake and in the intermediate wake, respectively, due to nonlinearities induced by the convective acceleration and steady/unsteady interactions. The nature of the unsteady displacement thickness suggests viscousinviscid interaction as a possible mechanism for the axial pressuregradient response.
Future work involves a parametric study for the ffx geometry for combined and verticalwave inflows (Appendix A) to determine the effects of frequency and extensions for realistic unsteady propeller/hull/appendage flows.
ACKNOWLEDGMENTS
This research was sponsored by ONR under Contracts N00014–92J1118 and N00014– 91J1203 under the administration of Mr. Jim Fein and Dr. Ed Rood whose support is greatly appreciated. Computer funds were provided by the NASA Numerical Aerodynamic Simulation Program and by ONR NAVOCEANO Primary Oceanographic Prediction System. The first author is grateful for the support provided by the IIHR Hunter Rouse Fellowship.
REFERENCES
1. Fage, A. and Johansen, F.C., “On the Flow of Air Behind an Inclined Flat Plate of Infinite Span, ” Proc. Royal Society A, Vol. 116, 1927, pp. 170–197.
2. Tyler, E., “Vortex Formation Behind Obstacles of Various Sections,” Philosophical Magazine, Vol. 11, No. 7, 1931, pp. 849–891.
3. Bearman, P.W. and Graham, J.M.R., “Vortex Shedding from Bluff Bodies in Oscillatory Flow: a Report on Euromech 119,” J. Fluid Mech., Vol. 99, pp. 225–245.
4. Shirayama, S., “Flow Past a Sphere: Topological Transitions of the Vorticity Field, ” AIAAJ., Vol. 30, No. 2, 1992, pp. 349–366.
5. Perry, A.E. and Watmuff, J.H., “The PhaseAveraged LargeScale Structures in ThreeDimensional Turbulent Wakes,” J. Fluid Mech., Vol. 103, pp. 38–51.
6. Theodorsen, T, “General Theory of Aerodynamic Instability and the Mechanism of Flutter, ” NACA Report No. 496, 1934.
7. McCroskey, W.J., “Unsteady Airfoils,” Ann. Rev. Fluid Mech., Vol. 14, 1982, pp. 285– 311.
8. Ghia, K.N., Yang, J., Ghia, U. and Osswald, G., “Analysis of Dynamic Stall Phenomenon through Simulation of Forced Oscillatory Flows,” Proc. 5th Symp. on Numerical and Physical Aspects of Aerodynamic Flows, Long Beach, CA, 1992.
9. Ohmi, K., Coutanceau, M., Loc, T.P., and Dulieu, A., “Vortex Formation around an Oscillating and Translating Airfoil at Large Incidences,” J. Fluid Mech., Vol. 211, 1990, pp.37–60.
10. Telionis, D.P., Unsteady Viscous Flow, New York, SpringerVerlag, 1981.
11. Commerford, G.L. and Carta, F.O., “Unsteady Response of a TwoDimensional Airfoil at High Reduced Frequency, ” AIAA J., Vol. 12, No. 1, 1974, pp. 43–48.
12. Griffin, O.M. and Hall, M.S., “ReviewVortex Shedding Lockon and Flow Control in Bluff Body Wakes, ” ASME J. Fluids Engineering, vol. 113, no. 4, 1991, pp. 526–537.
13. Katz, Y., Nishri, B., and Wygnanski, I., “The Delay of Turbulent Boundary Layer Separation by Oscillatory Active Control,” AIAA 2nd Shear Flow Conference (AIAA 89–0975), Tempe, AZ, 1989.
14. Favier, D., Maresca, Ch., and Castex, A., “Experimental Study of 2D/3D Interactions Between a Vortical Flow and a Lifting Surface,” Eur. J. Mech., B/Fluids, Vol. 8, No. 5, 1989, pp. 397–430.
15. Srinivasan, G.R., McCroskey, W.J. and Baeder, J.D., “Aerodynamics of TwoDimensional BladeVortex Interaction,” AIAA J., Vol. 24, No. 10, 1986, pp1569–1576.
16. Rai, M.M., “ThreeDimensional NavierStokes Simulations of Turbine RotorStator Interaction; Part IIResults,” J. Propulsion, Vol. 5, No. 3, 1989, pp. 312–319.
17. Huerre, P. and Monkewitz, P.A., “Local and Global Instabilities in Spatially Developing Flows,” Ann. Rev. Fluid Mech., Vol. 22, 1990, pp. 473–537.
18. Stern, F., Kim, H.T., Zhang, D.H., Toda, Y., Kerwin, J., and Jessup, S., “Computation of Viscous Flow around PropellerBody Configurations: Series 60 C_{B}=.6 Ship Model,” to appear, J. Ship Research.
19. Kim, H.T. and Stern, F., “Viscous Flow around a PropellerShaft Configuration with InfinitePitch Rectangular Blades,” J. Propulsion Vol. 6, No. 4, 1990, pp. 434–444.
20. Stern, F., Zhang, D.H., Chen, B., Kim, H.T., and S.D.Jessup, “Computation of Viscous Flow around PropellerShaft Configurations, ” to appear, Proc. 20th ONR Symposium on Naval Hydrodynamics, Santa Barbara, CA, August 1994.
21. Kim, W.J., Chen, B. and Stern, F., “Efficient Unsteady ThreeDimensional ViscousFlow Method with Application to Marine Propulsors,” to appear, Proc. 20th ONR Symposium on Naval Hydrodynamics, Santa Barbara, CA, August 1994.
22. Paterson, E.G., “Computation of Natural and Forced Unsteady Viscous Flow,” Ph.D. thesis, The Univerity of Iowa, expected May 1994.
23. Chen, H.C. and Patel, V.C., “Calculation of TrailingEdge, Stern and Wake Flows by a TimeMarching Solution of the Partially
Parabolic Equations,” IIHR Report No. 285, Iowa Institute of Hydraulic Research , 1985.
24. Rodi, W. and Srinivas, K., “Computation of Flow and Losses in Transonic Turbine Cascades,” Z.Flugwiss. Weltraumforsch, vol. 13, 1989, pp. 101–119.
25. Cebeci, T., Clark, R.W., Chang, K.C., and Halsey, N.B., “Airfoils with Separation and the Resulting Wakes,” Proc. 3rd Symp. on Numerical and Physical Aspects of Aerodynamic Flows, Long Beach, CA, 1985.
26. Granville, P.S., “BaldwinLomax Factors for Turbulent BoundaryLayers in Pressure Gradients, ” AIAA J., Vol. 25, No. 12, 1987, pp. 1624–1627 .
27. Roache, P.J., Computational fluid dynamics, Albuquerque, Hermosa publishers, 1972, pp. 184–185.
28. Thompson, J.F., “A Composite Grid Generation Code for General 3D Regions,” AIAA 25th Aerospace Sciences Meeting, Reno, NV, 1987.
29. Benek, J.A., “Extended Chimera Grid Embedding Scheme with Applications to Viscous Flows,” AIAA 8th Computational Fluid Dynamics Conference, 1987, Honolulu.
30. Suhs, N.E., and Tramel, R.W., “PEGSUS 4.0 User's Manual,” Arnold Engineering Development Center, AEDCTR91–8, 1991.
31. Rice, J.Q., “Investigation of a TwoDimensional Hydrofoil in Steady and Unsteady Flows,” M.S.Thesis, Massachusetts Institute of Technology, 1991.
32. Delpero, P.M., “Investigation of Flows around a TwoDimensional Hydrofoil Subject to a High Reduced Frequency Gust Loading,” M.S. Thesis, Massachusetts Institute of Technology, 1992.
33. Horwich Lurie, B. “Unsteady Response of a TwoDimensional Hydrofoil Subject to High Reduced Rrequency Gust Loading,” M.S.Thesis, Massachusetts Institute of Technology, 1993.
34. Chung, K.N., “Applications of a Zonal Approach and Turbulence Models to Foils and Marine Propellers,” Ph.D Thesis, The University of Iowa, December 1993.
35. Fuhs, D., “Comparisons of TwoDimensional Unsteady RaNS Calculations and Measurements, ” Carderock Divsison, Naval Surface Warfare Center report (in preparation) .
36. Taylor, L.K., Busby, J.A., Jiang, M.Y., Arabshahi, A., Sreenivas, K., Whitfield, D.L., “TimeAccurate Incompressible NavierStokes Simulation of the Flapping Foil Experiment,” Sixth International Conference on Numerical Ship Hydrodynamics, Iowa City, Iowa, Augst 1993.
37. Fan, S., Lakshminarayana, B., Barnett, M., “LowReynoldsNumber kƐ Model for Unsteady Turbulent BoundaryLayer Flows,” AIAA J., Vol. 31, No. 10, 1993, pp. 1777–1784.
38. Patel, M.H., “On Laminar Boundary Layers in Oscillatory Flow,” Proc. R. Soc. Lon. A. Vol. 347, 1975, pp. 99–123.
39. Choi, J.E., “Role of FreeSurface Boundary Conditions and Nonlinearities in Wave Boundary Layer and Wake Interaction,” Ph.D Thesis, The University of Iowa, December 1993.
Appendix: A:
Tunnel domain inflow, potential flow model
The inflow is divided into three regions: (1) between the bottom tunnel wall and the bottom flapper wake; (2) between the flapper wakes; and (3) between the top flapper wake and the top tunnel wall. For (1) and (3) an analogy can be made between the sinuous wake sheet of the oscillating flappers passing over the tunnel walls and a shallowwater wave. The potential flow describing these regions is:
The region between the flappers (although not used in this paper) is
A study was conducted to compare this model to two other models, one using unsteady, sinusoidal vortex sheets, and one using discrete Oseen vortices convecting at the freestream velocity. The chosen model was easiest to apply and gave the best correlation to the inlet data.
DISCUSSION
by Dr. Ming Zhu, University of Tokyo
First, I would like to congratulate your nice work on the unsteady NS simulation. Please show your Fourier analyses results of field data about the vortexfoil interaction. I would like to know that with your formulation, how much the unsteadiness can be caught.
Authors' Reply
The first twenty harmonics where calculated for the boundarylayer profiles which showed that the frequency content was dominated by the first harmonic. However, the wake of foil showed a much richer frequency content. As far as the temporal resolution of the method is concerned, it is only limited by the time step used in the simulation. Since the method is implicit, fairly large time steps are used, thus making the resolution of higher harmonics difficult. However, if a broadband or higher frequency response was expected, the time step could be reduced such that the frequencies of interest were sufficiently resolved.
DISCUSSION
by Dr. D.H.Choi, Korea Advanced Institute of Science and Technology
The calculations have been performed with three different computational domains. In principle, if the boundary conditions are correct, the small domain solution should be most accurate as it eliminates any uncertainty or inaccuracy associated with outer regions. However, since available experimental data are sparsely distributed, the boundary conditions for the small domain calculation are supplemented by using the tunnel domain calculation. This essentially makes these two calculations identical as the small domain is a subdomain of the tunnel domain. Therefore, the comparison of the two appears moot.
Authors' Reply
This is exactly the point that we wanted to make. At the ffx workshop, the participants who used the small domain had similarly poor solutions, i.e., erroneous higher harmonics in the pressure field. It was not clear that the small domain could be used for simulation of the ffx. Therefore, to demonstrate the validity of the smalldomain boundaryconditions, the tunnel domain solution was used in place of the data with the intention of making this moot comparison the proof that the formulation was correct.
DISCUSSION
by Dr. Charles C.Song, University of Minnesota
Did you compare the calculated unsteady pressure with experimental value? One should expect large differences between numerical results based on incompressible flow theory and the experiment, especially for higher frequency components.
Authors' Reply
Yes. The differences were large, but were most likely due to difficulties in obtaining accurate unsteady pressure data, not compressibility effects.
Time Accurate Incompressible NavierStokes Simulation of the Flapping Foil Experiment
L.K.Taylor, J.A.Busby, M.Y.Jiang, A.Arabshahi, K.Sreenivas, and D.L.Whitfield (CFD Laboratory, USA)
ABSTRACT
A numerical simulation of an experiment conducted in the Marine Hydrodynamics Laboratory at Massachusetts Institute of Technology (MIT) has been performed. The experiment was designed to study the flow about a twodimensional hydrofoil undergoing high reduced frequency gust loading. The gust was created by two NACA 0025 hydrofoils oscillating sinusoidally in phase upstream of the stationary foil. Experimental data was taken in the flowfield near the stationary foil as well as on the foil itself. The numerical computations were performed on the entire experimental domain, including the flapping foils, through the use of a twodimensional multiblock unsteady incompressible NavierStokes algorithm based on artificial compressibility. A comparison of the pressures and velocities on the bounding box and surface of the stationary foil are presented for the steady case, as well as time histories and an harmonic analysis at the same locations for the unsteady case.
NOMENCLATURE
A,B, 
Flux Jacobians 
Roe Matrix 

c_{p} 
Pressure coefficient 
E,F,K 
Computational space flux vectors 
ê, 
Nu^{m}erical flux vectors 
J 
Metric Jacobian 
m 
Newton iteration parameter 
n 
Time level index 
p 
Pressure 
Q 
Computational space dependent variable vector 
q 
Physical space dependent variable vector 
Re_{L} 
Reynolds number based on reference length 
T,T^{−1} 
Eigenmatrices 
t 
Time 
u,v 
Cartesian velocity components 
U,V,θ_{k} 
Contravariant velocities 
x,y 
Cartesian coordinates 
β 
Artificial compressibility parameter 
∂ 
Partial differentiation 
δ 
Central difference operator 
μ 
Viscosity 
ϱ 
Density 
τ 
Time in computational space; fluid stress vector 
ξ,η,k 
Curvilinear coordinates 
INTRODUCTION
Steady and unsteady computations were made of the MIT/ONR Flapping Foil Experiment (1,2) using the Mississippi State University incompressible NavierStokes code (3). This code is threedimensional; however, a twodimensional version was developed for these computations. The artificial compressibility form of the equations is solved on a timedependent curvilinear coordinate system. The equations are discretized in finite volume form and the numerical flux formulation at cell faces is based on Roe's approximate Riemann solver (4). The equations are formulated for numerical solution as a formal Newton method where the time derivative of the dependent variable vector, except for the time derivative of pressure in the continuity equation, is included in the residual term. This system is then solved at each time step using what is referred to as discretized Newtonrelaxation (5). Relaxation is carried out at each time step using three symmetric GaussSeidel passes. The computations are secondorder accurate in time,. For steady state solutions, local time stepping is used and the Jacobian matrix is
* Research Engineer ** Graduate Student + Professor 
updated every 20 cycles. A BaldwinLomax turbulence model is used for all the computations.
The NavierStokes code is a multiblock dynamic grid code and the computations are performed using the rather complete configuration consisting of the water tunnel walls, flapping foils, and stationary foil (Figs. 1 and 2). No slip boundary conditions are used on the tunnel walls, flapping foils and stationary foil. A steady state solution is first obtained and then motion of the flapping foils is initiated by boundary conforming dynamic grids that move in pitch with the flapping foils at a reduced frequency of 3.62 based on the halfchord of the stationary foil. The grid extends approximately two stationary foil chord lengths upstream of the leading edge of the stationary foil and five stationary chord lengths downstream of the trailing edge of the stationary foil. (Note that the entire grid is not presented in Fig. 2.) Characteristic variable boundary
conditions are used at both the upstream inflow boundary and downstream outflow boundary. Velocities are specified at the inflow while static pressure is set at the outflow.
For these multiblock computations, four blocks are used. Block one extends from the lower tunnel wall to the centerline of the lower flapping foil. Block two extends from the top of block one to the centerline of the stationary foil. Block three extends from the top of block two to the centerline of the upper flapping foil. Block four extends from the top of block three to the upper tunnel wall. All grids are Htype grids. This blocking arrangement allows the multiblock solutions to correspond exactly to a single block solution, even for unsteady flow. That is, there are no approximations at block boundary interfaces.
A presentation of the equations and numerical discretization will be given first, followed by a discussion of the numerical flux formulation and numerical solution method Then a section containing a description of the stationary and dynamic grids that are used followed by an explanation of the post processing procedures that are necessary to examine the results. Finally, results for both the steady and unsteady cases along with conclusions are presented.
EQUATIONS AND NUMERICAL DISCRETIZATION
Governing Equations
The unsteady twodimensional incompressible NavierStokes equations without body forces are written in Cartesian coordinates and in conservative form. The quantities used to nondimensionalize these equations are, a characteristic length, and freestream values of velocity, density, and viscosity. Pressure is normalized using the following relationship , where the subscript denotes a freestream or reference value (3). The resulting set of nondimensional equations are then transformed to a general timedependent bodyconforming curvilinear coordinate system given by
ξ=ξ(x,y,t)
η=η(x,y,t)
τ=t.
(1)
The timedependent nature of this transformation will allow all computations to be carried out on a fixed uniform computational domain even though components of the physical domain may be in motion.
Artificial compressibility is introduced after the transformation to insure that the continuity equation is satisfied for dynamic grid applications (6). The resulting set of equations may be written as
(2)
where
and
U=ξ_{t}+ξ_{x}u+ξ_{y}ν
V=η_{t}+η_{x}u+η_{y}ν
with
where Re_{L} is the Reynolds number based on the characteristic length. J is the Jacobian of the inverse transformation , given by
J=x_{ξ}y_{η}−y_{ξ}x_{η}
and the metric quantities are
Numerical Discretization
A finitevolume discretization of Eq. (2) in twodimensional computational space may be written as
(3)
where the indices i, j correspond to a cell center and the central difference operator δ is given by
. (4)
The increments Δξ and Δη are set equal to one for simplicity and Eq. (3) becomes
. (5)
In this expression, the dependent variable vector Q is considered to be constant throughout cell i, j, while the inviscid and viscous fluxes are assumed to be uniform over each of the four surfaces of the cell.
Temporal Discretization
The time derivative appearing in Eq. (5) is approximated with a general difference expression (7), which has been written for computational space, and is given by
(6)
where
ΔQ^{n}=Q^{n+1}−Q^{n}.
The two choices for θ_{1} and θ_{2} that were used exclusively in this work are θ_{1}=1, θ_{2}=0 and θ_{1}=1, θ_{2}=0.5, which correspond respectively to the backward Euler implicit scheme (firstorder accurate in time) and the three point backward implicit
scheme (secondorder accurate in time). Using Eq. (5) in Eq. (6) yields
(7)
where
R*=[δ_{i}(E−E_{v})+δ_{j}(F−F_{v})]^{*}.
In this particular application, the movement of the flapping foils dictate that a new grid must be generated at each time step. This implies that the Jacobian of the inverse transformation will be a function of the nondimensional time τ. The approach suggested by Janus (8) will be used to account for this dependence. First of all, the following two identities are needed:
ΔQ^{n}=Δ(J_{q})^{n}≡J^{n+1}Δq^{n}+q^{n}ΔJ^{n}
(8)
ΔQ^{n−1}=Δ(J_{q})^{n−1}≡J^{n−1}Δq^{n−1}+qnΔJ^{n−1}.
Inserting Eq. (8) into Eq. (7) and simplifying leads to
(9)
where
.
According to Janus (8). the evaluation of the bracketed term in Eq. (9) must be done with the same numerical scheme that is used to advance the NavierStokes equations in time in order to satisfy the geometric conservation law (9). Toward this end, consider the twodimensional discretized geometric conservation law,
(10)
Using the temporal discretization given by Eq. (6) yields
(11a)
where
.
(11b)
Now, the left hand side of Eq. (11a) is identically the bracketed term of Eq. (9), thus Eq. (9) can be written as
(12)
The evaluation of both the convective and diffusive flux vectors contained in R will be discussed in the next section.
NUMERICAL FLUX FORMULATION
Convective Flux Evaluation
FluxDifference Splitting
In this cellcentered finitevolume scheme, the convective numerical flux at cell faces is calculated using the Roe approximate Riemann solver (4). Although originally developed for compressible flow, this approximate Riemann solver can also be implemented in artificial compressibility formulations for incompressible flow (additional details may be found in Reference (10)). An essential ingredient of Roe's solver is the construction of a matrix , which is representative of the local interface conditions. Roe (4) required this matrix to satisfy certain specific properties, two of which are
.
The satisfaction of all these properties in compressible flow is guaranteed by uniquely averaging (11) the dependent variables and then using them in the evaluation of the matrix These specially averaged variables are commonly referred to as Roe variables. An analogous set of Roe variables can be defined for incompressible flow (6,12,13) and are given by
Once these variables have been defined, the mechanics involved in computing the numerical fluxes are the same as those used for compressible flow.
A formulation of the numerical flux at cell face i+1/2 which uses all Roe variables and metrics
corresponding to cell face i+1/2 can be found in Reference (14). This higherorder numerical flux without limiters is
(13)
where
and
λ^{±}corresponds to the positive and negative eigenvalues of the Roe matrix and r^{(j)} and ℓ^{(j)} are the right and left eigenvectors, respectively. The subscript i+1/2 in the above equations indicate that the metrics used are evaluated at cell face i+1/2. The dependent variables in the eigenvalues and eigenvectors are the aforementioned Roe variables. The flux vector e(q_{i}), however, is evaluated using the dependent variables. All of the results reported herein were obtained using this numerical flux vector with ψ=1/3, which is thirdorder accurate in space.
Eigensystem of the Transformed Flux Jacobian Matrices
It is apparent from Eq. (13) that the eigenvalues and left and right eigenvectors of the Roe matrix must be known in order to perform the numerical flux calculation. The eigensystem for the Roe matrix can be obtained from the eigensystem of the transformed flux Jacobian by simply using the Roe variables in the latter.
The first step in deriving this eigensystem is the determination of the transformed flux Jacobian matrices. These matrices are defined as
(14)
where E and F are given in Eq. (2) and the vector q is
.(15)
The inviscid flux vectors E and F can be written compactly as
(16)
where
θ_{k}=k_{t}+k_{x}u+k_{y}ν
and
K=E, θ_{k}=U for k=ξ
K=F, θ_{k}=V for k=η.
Expressing K in terms of the elements of q and performing the indicated differentiation as defined in Eq. (14) yields the transformed flux Jacobian matrices,
(17)
where
for m=x,y, and t
and
.
To facilitate the determination of the eigenvalues of the matrix , a similarity transformation will be introduced, which will operate on and produce the matrix x. Since and x are similar matrices, they will have the same characteristic polynomial (15) and hence the same eigenvalues. However, the matrix x will contain fewer nonzero elements than matrix ,
and thus the computation of its eigenvalues will be far easier.
Consider the matrices M and M^{−1}
(18) .
Let , then
(19)
where
.
The eigenvalues of x can be calculated by using the first row to expand the determinant, which yields
(20)
where
and , i=1,2,3 are the eigenvalues of A for k= ξ and of B for k=η.
The right eigenvectors of x are solutions corresponding to the respective eigenvalues of the equation
, s=1,2,3.
(21)
The matrix P_{k,} whose columns are the right eigenvectors of x, can be written as
(22)
where
Having developed a linearly independent set of right eigenvectors, the left eigenvectors will be constructed in such a way that the right and left eigenvectors will be orthonormal. One way to accomplish this task is to simply invert the matrix P_{k}, which yields
(23) .
The columns of P_{k} and the rows of are, respectively, the right and left eigenvectors of x corresponding to the particular eigenvalues. In fact x can be written as
(24)
where Λ_{k} is a diagonal matrix whose elements are the eigenvalues λ_{k.} Recall that the original intent was to develop the eigensystem of the flux Jacobian matrices A and B. Having generated the eigensystem of x, the eigensystem of the flux Jacobians can be obtained as follows. Using Eq. (19), can be expressed as
.
(25)
Substituting the expression for x in Eq. (24) into Eq. (25) yields
(26)
or
(27)
where
Performing the indicated multiplication and simplifying, results in the following expressions for T_{k} and
(28)
Viscous Flux Evaluation
In this work, all velocity derivatives needed in the evaluation of the viscous flux vectors E_{v} and F_{v} are approximated using central differences. Cross derivative terms are calculated by simply averaging the appropriate metric quantities and dependent variables from surrounding cells. Gatlin (16) has given a detailed account of the evaluation of the remaining terms and it will not be repeated here.
NUMERICAL SOLUTION METHOD
The numerical solution of the nonlinear system of equations given by Eq. (12) is obtained by using Newton's method. This formulation is used for both steady state and unsteady solutions. The iterative nature of Newton's method will insure the compatibility of the computed pressure field and the divergencefree velocity field for unsteady calculations. The Jacobian matrix used in this solution procedure is computed by numerically differentiating the firstorder Roe flux vector. This discretized Jacobian is then used in a Newtonrelaxation scheme where Newton is the primary iteration and relaxation is the secondary iteration. Formally, Ortega and Rheinboldt (5) refer to this particular solution methodology as a discretized Newtonrelaxation method, or DNR (17). Reference (18) contains the basic derivation of the aforementioned scheme, while Reference (10) describes the implementation for incompressible flow. The remainder of this section will therefore only touch upon some of the fundamental concepts, while providing the details of its implementation for this particular application.
Newton Formulation and Jacobian Evaluation
Consider a system of nonlinear equations given by
N_{1} (x_{1},x_{2}, . . . , x_{n})=0
N_{2} (x_{1},x_{2,} . . . , x_{n})=0
.
.
.
N_{n} (x_{1,}x_{2}, . . . , x_{n})=0 (29)
or more concisely
N(x)=0. (30)
Newton's method for the vectorvalued function N(x) can be expressed as
N'(x^{m})=(x^{m}^{+1}−x^{m})=−N (x^{m}) (31)
where m=1, 2, 3, … and N′ (x^{m}) is the Jacobian matrix of N(x^{m}) given by
(32)
where . The approach used to solve Eq. (31) is to replace the analytical Jacobian with a numerically approximated Jacobian. Ortega and Rheinboldt (5) refer to such a method as a discretized Newton iteration. Various finitedifference approximations have been suggested in the numerical analysis literature. A simple and straightforward approximation of the elements of the Jacobian is to replace the analytical Jacobian with
where e_{j} is the j^{th} unit vector. Dennis and Schnabel (19) point out that this is the same as approximating the j^{th} 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 q^{n+1}, that is, in Eq. (31) x=q^{n+1}. Thus, Eq. (31) becomes
N′ (q^{n+1,m})(q^{n+1,m+1}−q^{n+1,m})=−N(q^{m}) (33)
where now the a_{ij} in Eq. (32) are given by
(34)
It should be noted that the only terms in Eq. (12) that are functions of q^{n+1} are Δq^{n} and R^{n+1}. Recall that R^{n+1} is defined as
R^{n+1}=[δi(E−E_{v})+δ_{j}(F−F_{v})]^{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 firstorder contribution to the higherorder 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 crossderivative 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 lefthand side of Eq. (33) is
(38)
where
and
Note that the first subscript of the Dg (g=e,e_{v},f or f_{v}) 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 righthand side of Eq. (33) is
(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 I_{a} 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 firstorder accurate, which leads to the definition of the matrices I_{b} and I_{c} 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 GaussSeidel 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 GaussSeidel sweeps per Newton iteration.
Turbulence Modeling
The algebraic turbulence model of BaldwinLomax (20), which was used in Reference (3) for threedimensional incompressible flow has been modified for twodimensional flow. The turbulent wake model used is that of Renze, Buning and Rajagopalan (21). Multiple peaks in the F_{max} 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 threedimensional grids has been incorporated into a twodimensional 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
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 blockstructured 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 blockstructured grid techniques. The blocked grid idea is not new. It has been applied to threedimensional 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 blockblock 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 blockblock interface communication is a direct extractioninjection 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 Htype 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
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 W_{i} and W_{j}, 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−W_{ij}) X(fix)+W_{ij}Lcos(α)
Y(new)=(1−W_{ij}) Y(fix)+W_{ij}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.
POSTPROCESSING
The computed results had to go through an elaborate postprocessing procedure in order to be analyzed in an accurate and meaningful manner. In
particular the computed results had to be rereferenced 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 rereference 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 steadystate 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 steadystate computations, the flapping foils were held fixed. Since a steady state solution
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 GaussSeidel 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 ucomponent 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.
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 uvelocity 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 BaldwinLomax 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.
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 semichord 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 uvelocity 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.
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 uvelocity 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.
CONCLUSIONS
A twodimensional multiblock unsteady incompressible NavierStokes 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 –92J1060 with James A. Fein as the technical monitor. The authors wish to express their sincere appreciation to the following organizations for providing CRAYYMP computer resources: U.S. Army Waterway Experiment Station and Cray Research, Incorporated.
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 ThreeDimensional Incompressible Euler and NavierStokes Solver for Stationary and Dynamic Grids,” AIAA91–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,” AIAA89–0122, January 1989.
7. Beam, R.M., and Warming, R.F., “An Implicit Factored Scheme for the Compressible NavierStokes Equations, ” AIAA Journal, Vol. 16, No. 4, pp. 393–402, April 1978.
8. Janus, J.M., “Advanced 3D 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 ThreeDimensional 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 NavierStokes Solutions for a SharpEdged DoubleDelta Wing,” AIAA 87–0206, January 1987.
13. Rogers, S.E. and Kwak, D., ”Upwind Differencing for the TimeAccurate Incompressible NavierStokes 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 WaveSplit Scheme for Solving the Unsteady ThreeDimensional Euler and NavierStokes Equations on Stationary or Dynamic Grids,” Engineering and Industrial Research Station Report, MSSUEIRSASE88 –3, Mississippi State University, Mississippi State, February 1988.
15. Sokolnikoff, I.S. and Redheffer, R.M., Mathematics of Physics and Modern Engineering, McGrawHill, Inc., New York, NY, 1966.
16. Gatlin, B., “An Implicit, Upwind Method for Obtaining Symbiotic Solutions to the ThinLayer NavierStokes Equations,” PhD Dissertation, Mississippi State University, August 1987.
17. Vanden, K.J., “Direct and Iterative Algorithms for the ThreeDimensional Euler Equations, ” PhD Dissertation, Mississippi State University. December 1992.
18. Whitfield, D.L., and Taylor, L.K., “Discretized NewtonRelaxation Solution of High Resolution FluxDifference Split Schemes,” AIAA91–1539, June 1991.
19. Dennis, J.E., Jr. and Schnabel, R.B., Numerical Methods for Unconstrained Optimization and Nonlinear Equations , PrenticeHall, Inc., Englewood Cliffs, New Jersey, 1983.
20. Baldwin, B.S. and Lomax, H., “ThinLayer Approximation and Algebraic Model for Separated Turbulent Flows,” AIAA78–257, January 1978.
21. Renze, K.J., Buning, P.G. and Rajagopalan, R.G.,” A Comparative Study of Turbulence Models for Overset Grids,” AIAA92–0437, January 1992.
22. Chen, J.P., “Unsteady ThreeDimensional ThinLayer NavierStokes 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, ” AIAA84–0164, AIAA 22nd Aerospaces Sciences Meeting, Reno, NV, January 1984.
25. Thompson, J.F., “A Composite Grid Generation Code for General ThreeDimensional Regions, ” AIAA87–0275, AIAA 25th Aerospace
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,” AIAA84–1665, June 1984.
27. Belk, D.M., “ThreeDimensional 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. AU17, 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 BaldwinLomax 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 pressuregradient modifications of the BaldwinLomax turbulence model have been included in the paper. It is important to emphasize that Granville's change in the zerogradient values of C_{kleb} and C_{cp} had larger influence in improving the solution than the pressuregradient corrections to these coefficients.