Engineering Analysis with Boundary Elements 24 (2000) 277–286 www.elsevier.com/locate/enganabound
A simple and effective zonal procedure for lifting surfaces T.J. Barber*, E. Leonardi, R.D. Archer University of New South Wales, Sydney 2052, Australia Received 22 March 1999; received in revised form 18 October 1999; accepted 17 November 1999
Abstract Accurate prediction of viscous flows is currently achieved with a Navier Stokes solver, however these are computationally demanding. Panel methods are very efficient in predicting inviscid flow characteristics, however they are less accurate in regions of highly viscous flow. A method combining a panel method with a Navier Stokes solver has been developed, however unlike other zonal procedures, the current method is a simple, non-iterative linking procedure. Significant time and memory savings were found when compared with the use of the CFX (Navier Stokes) solver alone. 䉷 2000 Elsevier Science Ltd. All rights reserved. Keywords: Navier Stokes solver; Inviscid flow; Zonal procedures
1. Introduction Following Prandtl [1], Chiu [2] confirms that for flow computations at high Reynolds numbers, such as that for flow over vehicles, the flow field can “be perceived as large domains of essentially inviscid flow and comparatively smaller regions in which viscous effects are important.” The current work has been investigated during a study of wing-in-ground effect vehicles, for which typical Reynolds numbers are around 10 million and for which the flow-field can be so described. To solve such a flow accurately, a Navier Stokes solver can be used. However, large numbers of grid points are required in order to accurately predict flow in high gradient areas and for the computational boundaries to extend far enough from the body. Wendt [3] describes such a typical grid: “the white speck in the middle of the figure is the airfoil, and the grid spreads far away from the airfoil in all directions. The freestream is subsonic, hence the outer boundary must be placed far away from the airfoil because of the far-reaching propagation of disturbances in a subsonic flow.” An example of a typical Navier Stokes solution grid is shown in Fig. 1. The method proposed here uses a grid with boundaries no more extreme than those shown in Fig. 2. The effect of this proposed grid system would be to decrease substantially the amount of time required for the * Corresponding author. Tel.: ⫹ 612-9385-4093; fax: ⫹ 612-96631222. E-mail address:
[email protected] (T.J. Barber).
solution, and the amount of computer resources required to solve and save the problem. The proposed procedure is described in Fig. 3. It is based on the assumption that the flow may be divided into regions of viscous and inviscid flow. The procedure can be explained as follows. First, the grid is developed, using one of the pre-processing packages of CFX [4]. Boundaries at the upper, lower and left sides of the grid are defined as Dirichlet boundaries (or “Inlet”), allowing a specification of inlet boundary conditions (in this case, velocities). The CFX Fortran user subroutines examine the grid and extract only the grid points which form the surface of the lifting body, as panel methods have the advantage of requiring only the surface geometry of a body to be analysed. This information is passed to the panel method, which is then used to solve the singularity distribution on the body surface (see Section 2.1). Once this is known, the velocity at any off-body point can be calculated, by considering the distance from, and the strength of each singularity, and the free-stream conditions. The grid is again examined, this time, grid points which form the upper, lower and left side boundaries (those specified as Inlet) are extracted. The panel method is then used to calculate the velocity at these boundary points due to the effect of the lifting body. These velocities are provided as inlet conditions to the CFX solver, which then uses the standard solution process (see Section 2.2) to iterate to a required level of convergence. Although the concept of combining two such methods has been used before, previous attempts have used a complex
0955-7997/00/$ - see front matter 䉷 2000 Elsevier Science Ltd. All rights reserved. PII: S0955-799 7(99)00053-3
278
T.J. Barber et al. / Engineering Analysis with Boundary Elements 24 (2000) 277–286
Nomenclature B c C Cd Cl Cp r S t u U V v a G m r s sS F
body force chord body boundary coefficient of drag coefficient of lift coefficient of pressure distance from a point to a singularity source or sink term time variance of velocity time averaged velocity time averaged scalar variable variance of scalar variable angle of attack (degrees) diffusion coefficient doublet strength density stress tensor source strength velocity potential
system of iteration in order to allow the procedure to be used for varied applications and for the absolute minimization of zone boundaries. The current procedure uses a simple, noniterative approach to link the two methods. For example, Cebeci et al. [5] combined a panel method with a boundary layer method, iterating between the two solutions. The procedure begins by assuming a displacement thickness and solves both the inviscid and viscous parts of the solution using this value. Values for the velocity from both methods are calculated, and a relaxation formula is introduced to define a new displacement thickness, iterating until the two solutions converge. The results show improvement on previous inviscid procedures. Similarly, Navier Stokes solvers have been linked with a potential flow code. Roberts [6] used a panel method code, PMARC, to calculate velocity vectors and total pressure and temperature at the zonal crossover boundaries. Normal
Fig. 2. Proposed grid.
velocity boundary conditions were extracted from the Navier Stokes solution for use in PMARC. Summa et al. [7] also combined a panel method with the Navier Stokes method. The two computational zones overlap and after the first iteration, the potential flow surface panel strengths are obtained from the Navier Stokes solution at a smooth inner fluid boundary. The authors suggest the grid domain size can be reduced to about 0.14 chord lengths, with a 1–2% loss in accuracy. In 1995, Tuncer et al. [8] also used an overlapped zonal method but partitioned the computational domain into two zones—near field and farfield. The computational domain could be limited to a region that extended one-fifth of a chord length from the airfoil surface. The method was found to increase computational efficiency (in terms of CPU time) by 40%. The present procedure, however, aims to use a simple interpretation of the concept, specialized for a particular application, in order to exploit fully the time and cost savings. The investigation benefiting from this procedure considers airfoils and wings and so a basic panel method is ideal for the inviscid part of the solution. However, a detailed study of these surfaces near the ground requires accurate knowledge of the viscous effects and so a Navier Stokes solution with a fine grid is needed. The present procedure enables the grid to be sufficiently fine where needed and still to be solved in reasonable time, using reasonable computer memory. No iteration, other than the CFX solver, is required. Although the main section of the paper deals with a twodimensional procedure for an airfoil, results are also given for a case where ground effect is present and for a threedimensional wing. 2. Two-dimensional studies with CFXZone
Fig. 1. Typical Navier Stokes grid.
In order to test the procedure, a method was first developed for calculating the flow around airfoils in free air. Validation against experimental results could then be performed. Once this showed the process to be satisfactory, the method was extended to include ground effect.
T.J. Barber et al. / Engineering Analysis with Boundary Elements 24 (2000) 277–286
279
Fig. 3. Flow-chart of procedure.
2.1. The two-dimensional panel method For flow fields, which can be considered incompressible and inviscid, the potential flow equations can be used to determine the flow characteristics. In particular, Laplace’s equation can be used to find solutions for the strengths of singularity distributions placed on the body surface and these, together with an appropriate boundary condition allow velocity distributions to be calculated. The panel method used in this work is a constant potential method, which uses locally constant strength sources and doublets, distributed over the body surface. The body surface is divided into a number of sections or “panels” with each having a local source and doublet strength. A Dirichlet boundary condition is enforced, by specifying a constant internal total potential for the body. The boundary condition for the procedure can be written as 1 Z 2 ⴱ F i
x; z s ln r ⫺ m
ln r dC ⫹ F ∞ 2p C S 2n constant
where F iⴱ
x; z represents the value for the internal potential of the body. Source values are represented by s S, while doublet values are represented by m . C represents the boundary of the body, and r defines the distance from a particular surface point to a particular singularity. The velocity potential due to the source and doublet distributions is found and velocity values within the flowfield are calculated by considering the contribution from the singularity distribution at a particular point. By knowing both the singularity strength and the induced velocity, these can be summed to produce a value for any point a distance r, from the body.
2.2. Navier Stokes solver The general-purpose computer code CFX has been used
1
Fig. 4. Grid one.
Fig. 5. Chordwise Cp distribution for a ⫺0:6⬚:
280
T.J. Barber et al. / Engineering Analysis with Boundary Elements 24 (2000) 277–286
Fig. 6. Chordwise Cp distribution for a 2:9⬚:
Fig. 8. Chordwise Cp distribution for a 12:4⬚:
to solve the Reynolds-Averaged Navier Stokes equations: 7·U 0
2
close the set of equations, and the flow is assumed to be incompressible. 2.3. Numerical examples and validation of procedure
2U ⫹ 7·
U 丢 U B ⫹ 7·
s ⫺ u 丢 u 2t
3
2V ⫹ 7·
U V 7·
G 7V ⫺ uv ⫹ S 2t
4
These three equations represent continuity, momentum and scalar transport, where U is the time averaged velocity, s is the stress tensor, V is a mean scalar variable, G is the diffusion coefficient, B is a body force, S is a source or sink term and t is time. The terms, ru 丢 u; ruv represent Reynolds stress and Reynolds flux. The set of equations was discretized using a finite volume method and the SIMPLEC algorithm was used for pressure– velocity coupling. The Van Leer higher order differencing scheme was used to discretize the equations and the Block Stone solver was employed to solve the complete set of equations. The RNG k-e turbulence model was used to
In order to find the appropriate position of the outer boundary locations for the model, a full sized Navier Stokes solution grid was run with CFX. A NACA 4412 airfoil section at 6.4⬚ angle of attack and a Reynolds number of 1,800,000 were chosen. The grid had boundaries located at 10c (10 chord lengths) downstream, 6c upstream and 8c above and below the airfoil. Once this solution had been found, the panel method part of the CFXZone procedure was used to calculate velocity values at certain points within the grid, for specified distances in front of, above and behind the airfoil surfaces. These calculated values were compared with the CFX results at the same locations, and the maximum error found. • At 1/4c from the airfoil surface, the error in u-velocity values was found to be 7% and the error in v-velocity was found to be 8%. • At 2c from the airfoil surface the error in u-velocity was found to be 0.1% and the error in v-velocity was found to be 0.2%. • At 4c from the airfoil surface, the error in u-velocity was found to be 0.03% and the error in v-velocity was found to be 0.05%. 2.4. Test of symmetry
Fig. 7. Chordwise Cp distribution for a 6:4⬚:
In order to test that the procedure retains symmetry, a case was run with a NACA 0012 airfoil (a symmetrical section) at 0⬚ angle of attack. The expected solution for a case such as this is for the flow characteristics to be identical on the upper and lower surfaces of the airfoil. As the panel method solutions are formed as a direct result of geometric values, it was important to ascertain if small differences in the grid locations would affect the solution. The grid for this
T.J. Barber et al. / Engineering Analysis with Boundary Elements 24 (2000) 277–286
281
Fig. 9. Farfield grids.
case was chosen to be asymmetric, utilizing a different number of grid points on the upper and lower surfaces of the airfoil. The boundaries were located at 2c in front, above and below the airfoil and 3c behind the airfoil. Comparisons
were made between the calculated pressure coefficient distributions for the upper and lower surfaces of the airfoil. Results were found to be the same within machine round off error.
282
T.J. Barber et al. / Engineering Analysis with Boundary Elements 24 (2000) 277–286
Fig. 10. Cp distribution, farfield variations.
2.5. NACA 4412 airfoil at 6.4⬚ angle of attack A further test case used to validate the procedure was a NACA 4412 airfoil, for a Reynolds number of 1,800,000, and angles of attack of ⫺0.6, 2.9, 6.4 and 12.4⬚. The results were compared with the experimental data of Pinkerton [9]. The grid used has 3200 cells (Fig. 4). The grid points are clustered around the leading edge of the airfoil, as this is where the highest gradients occur. One hundred points form
the airfoil surface. The farfield boundary is situated 1/2c in front of the airfoil, and 1c above, below and behind the airfoil. Clearly, this is an extremely confined grid, and unsuitable for a standard Navier Stokes solution. Results are given for each case in the form of Cp curves, comparing the accuracy of using the current combined procedure, here called CFXZone, against CFX alone for the same grid. Whenever CFX is used, inlet conditions are specified as the free-stream flow conditions.
T.J. Barber et al. / Engineering Analysis with Boundary Elements 24 (2000) 277–286
283
Fig. 11. Comparison of CFX versus CFXZone, Cl, Cd.
For the a ⫺0:6⬚ case, CFXZone gives a better result than CFX for the upper and lower surfaces of the airfoil (Fig. 5). In all regions, the CFXZone result is very close to the experimental value. The CFX results, although predicting the correct trend, are not accurate. The inaccuracy in prediction can be seen most notably at the leading edge suction peak, where the experimental value is C p ⫺0:47 and CFX predicts a value of ⫺0.76. For a 2:9⬚; CFXZone is clearly more accurate than CFX (Fig. 6). The CFXZone results are accurate for all regions, while CFX only predicts the trailing edge values with accuracy. Leading edge values are notably inaccurate. For a 6:4⬚; CFXZone is again clearly more accurate than CFX (Fig. 7). Both pressure peaks are predicted well. CFX is inaccurate for most of the airfoil region, predicting the suction peak very poorly. For a 12:4⬚; CFXZone results show some inaccuracies for the latter part of the upper surface of the airfoil, however both pressure peaks are predicted well (Fig. 8). The reason for these inaccuracies is the use of the RNG k-e turbulence model, which is well known for its inability to cope with separated flow regions. CFX is inaccurate for
most of the airfoil region, predicting the suction peak very poorly. 2.6. Farfield boundary variation In order to determine the optimum position for the farfield boundaries when using CFXZone, various sizes of computational domain were used to calculate the pressure distribution for the a 6:4⬚ case. Five grids in total were used, with the farfield boundaries at the distances of 1.5c, 2c, 3c, 5c and 10c (cases 1–5). The grids are shown in Fig. 9. A CFX analysis was performed on each grid, and a CFXZone analysis was also performed on each of the grids. Results for the investigation are shown in Fig. 10 in the form of pressure coefficient distributions. Experimental results are also given for comparison [9,10]. The pressure distributions clearly show that the CFX results become successively more accurate as the boundaries are moved further away. The results calculated by CFX on the first four grids are however quite inaccurate. It is not until the fifth grid (with the boundaries located at 10c from the airfoil surface) is
Fig. 12. Efficiency comparisons.
284
T.J. Barber et al. / Engineering Analysis with Boundary Elements 24 (2000) 277–286
Fig. 13. Chordwise Cp distribution for ground effect case.
used that the results are accurate, compared to the experimental values. The results calculated by CFXZone however, show no such inaccuracies. Even for the first grid, the results compare very well with the experimental results. Further movement of the boundaries has little effect on the resulting pressure distribution. However, a comparison of the lift and drag values shows that these parameters do have some dependence on the boundary location. Comparing the values for CFX and CFXZone as the boundaries are moved further away (Fig. 11) it can be seen that CFX appears to require some further boundary distance to fully reach the converged value of lift and drag coefficients. CFXZone however shows convergence for these parameters by the third case, with the boundaries located at 3c from the airfoil surface. Comparing the computing workspace, CPU time, and iterations required for convergence shows the efficiency of CFXZone. It should also be noted that for the last three CFX grids some under-relaxation was required for CFX to reach convergence. No under-relaxation was required for any of the CFXZone cases. Comparing CFXZone (grid 3) and CFX (grid 5), (where the CFX solution is still not as accurate as the CFXZone solution), the total workspace has increased by ⬃30% for the CFX result. CPU time has increased by ⬃150%, and the total number of CFX iterations has increased by ⬃50% (Fig. 12). 3. NACA 4412 airfoil in ground effect In order to use the procedure for the application of ground effect, it was necessary to incorporate an added feature into the calculation. In the panel method section of the procedure this was achieved by using the “image” method, in which the panel
method reads in the surface points and then produces an identical set of points, except for having the y-coordinate negative. This gives the effect of calculating the flow over two identical airfoils, with a plane of symmetry at y zero. This plane of symmetry represents the ground. This is by far the easiest boundary condition to implement in a panel method and is an adequate representation for this application. Results were validated by comparing with a CFX solution calculated on a standard full-sized grid, with farfield boundaries positioned 10–20 chord lengths from the body. Results for CFXZone are shown for the NACA 4412 airfoil flying at a Reynolds number of 1,800,000 and an angle of attack of 6.4⬚. Ground clearance is 5% of the chord (Fig. 13). Clearly, the Cp curves are almost identical, showing that even for this extremely close case of ground effect CFXZone gives excellent results.
4. NACA 0012, AR 6 wing It would be especially advantageous to be able to limit the number of grid points required for a three-dimensional solution and therefore the procedure has been extended to three dimensions. A three-dimensional panel method code, based on the work of Hess and Smith [11], was written for this part of the work. A procedure similar to that of the two-dimensional method was used; the panel method reads a file of surface points, calculates a source and doublet distribution for the body, and then writes a file of velocity values for the grid boundary points. The test-case consists of a semi-span NACA0012 wing, aspect ratio 6, flying at angles of attack of 2.53, 6.75 and 13.06⬚ and at a Reynolds number of 3,500,000. Representative results for this test-case are shown in Fig. 14a–f. The results compare Cp curves for six cases: each of the three angles of attack at two spanwise locations. The results show excellent agreement with the experimental results [12] and the trend of decreasing suction peak values as the wing tip is approached is well represented.
5. Conclusions A simple aerodynamic analysis procedure has been developed, combining a commercial Navier Stokes solver (CFX) with a panel method code. The procedure has been developed for the particular case of airfoils and in modifications to the procedure, for airfoils with ground effect and for wings. The procedure, CFXZone, has been shown to be very cost effective, in both time and computer memory savings. The two-dimensional case was found to produce CPU time savings of 150%, workspace savings of over 30% and the
T.J. Barber et al. / Engineering Analysis with Boundary Elements 24 (2000) 277–286
Fig. 14. Chordwise Cp distribution comparison for three-dimensional procedure.
285
286
T.J. Barber et al. / Engineering Analysis with Boundary Elements 24 (2000) 277–286
number of iterations required for convergence was decreased by 50%. An extension to the procedure to allow ground effect calculations was found to provide excellent accuracy when compared to a full-sized grid CFX solution. Similarly, a three-dimensional extension of the new method showed excellent agreement with experimental results. References [1] Prandtl L. On fluid motions with very small friction. Proc Third Int Math Cong Heidelberg, 1904. English translation, NACA TM 452. [2] Chiu T, Broers C, Wood B. The application of modern panel method techniques to sport aero/hydrodynamics. AIAA Paper 95-2210, 1995. [3] Wendt J. Computational fluid dynamics. 2nd ed. New York: Springer, 1996.
[4] CFX-4.2. CFX International. AEA Technology. [5] Cebeci T, Sedlock D, Chang K, Clark R. Analysis of wings with flow separation. J Aircraft 1989;26(3):214–20. [6] Roberts D. A zonal method for modeling powered-lift aircraft flow fields. NASA CR 177521, 1989. [7] Summa J, Strash D, Yoo S. Zonal flow analysis method for twodimensional airfoils. AIAA J 1992;30(2):548–9. [8] Tuncer I, Ekaterinaris J, Platzer M. A novel viscous–inviscid interaction method with inviscid wake modeling. AIAA Paper 94-0178, 1994. [9] Pinkerton R. The variation with Reynolds number of pressure distribution over an airfoil surface. NACA Report No. 613, 1938. [10] Abbott I, von Doenhoff A. Theory of wing sections. New York: Dover, 1957. [11] Hess J, Smith A. Calculation of potential flow about arbitrary bodies. Prog Aero Sci 1966;8:1–138. [12] Yip L, Shubert G. Pressure distribution on a 1- by 3-metre semispan wing at sweep angles from 0⬚ to 40⬚ in subsonic flow. NACA TN D8307.