AERODYNAMIC LOADS OVER ARBITRARY BODIES BY METHOD OF INTEGRATED CIRCULATION by Vivek Ahuja A dissertation submitted to the Graduate Faculty of Auburn University in partial fulfillment of the requirements for the Degree of Doctor of Philosophy Auburn, Alabama August 3, 2013 Keywords: Potential-flow, Unstructured mesh, aerodynamic loads, Vorticity based drag prediction, circulation based loads, aerodynamic optimization Copyright 2013 by Vivek Ahuja Approved by Roy Hartfield, Chair, Woltosz Professor, Department of Aerospace Engineering John Burkhalter, Co-chair, Professor emeritus, Department of Aerospace Engineering Andrew Shelton, Assistant Professor, Department of Aerospace Engineering ii ABSTRACT A method for the application of vorticity based potential-flow solvers to unstructured surface meshes has been created. This method is designed to maintain the advantages of the vorticity based solvers while removing the limitations of geometrical inputs associated with the philosophical approach. A discussion on the necessity and advantages of this approach has been presented. A modification to the evaluation of skin- friction coefficients using surface vorticity has also been developed. An unstructured wake-strand model has been developed to allow handling of wakes emanating from unstructured meshes. An attempt has been made to extend potential-flow solvers to the current industry surface meshing and numerical solver standards. This automated pipeline is considered to be robust enough to allow integration with optimizers unrestricted in their geometrical design spaces. A test-suite of basic shapes and bodies was tested to evaluate the fidelity of the new approach. An advanced test-suite was developed to stress-test the solver as well as to test out the geometry handling pipeline required to handle these test cases. It was determined that the solver is able to provide high fidelity results for a wide variety of test cases and is able to: work with conformal and non-conformal geometry interfaces, resolve surface intersections, work with both structured and unstructured meshes, work with both thick and thin bodies and work with manifold and non-manifold surface interfaces. iii ACKNOWLEDGMENTS The author would like to thank Dr. Roy Hartfield for providing the encouragement, support and a free-thinking environment for conducting this research. Thanks are also due to Dr. John Burkhalter and Dr. Andrew Shelton for helping develop the fundamental theories used in this work and to Dr. Mark Carpenter for providing an objective overview of the work as an external reader of this document. Thanks are also due to my mother, Dr. Sulekha Ahuja, for providing that metaphorical wall of support to rest on during these past few wonderful years of work. iv Style manual or journal used: American Institute of Aeronautics and Astronautics Journal Computer software used: TECPLOT 10/360, FORTRAN v TABLE OF CONTENTS ABSTRACT .............................................................................................................................. ii ACKNOWLEDGMENTS ....................................................................................................... iii STYLE MANUAL................................................................................................................... iv LIST OF ILLUSTRATIONS .................................................................................................. vii LIST OF TABLES ................................................................................................................. xiv NOMENCLATURE ................................................................................................................xv 1. INTRODUCTION ...............................................................................................................1 2. MOTIVATION ................................................................................................................. 13 3. THE NEED FOR UNSTRUCTURED MESHES .............................................................16 4. SOLVER METHODOLOGY ............................................................................................26 4.1 Spatial velocity inductions ...........................................................................................26 4.2 Solver philosophy ........................................................................................................30 4.3 Vorticity wakes ............................................................................................................33 4.4 Aerodynamic loads ......................................................................................................41 4.4.1 Method of Integrated Circulation ..................................................................41 4.4.2 Skin friction drag evaluations .......................................................................50 5. MESHER METHODOLOGY ...........................................................................................53 5.1 Size based meshing ......................................................................................................53 5.2 Mesh and facet quality .................................................................................................59 5.3 Mesh auto-refinement ..................................................................................................61 5.4 Patch interface and conformity ....................................................................................69 vi 6. SOLVER EFFICIENCY ENHANCEMENT ...................................................................73 6.1 Fluid-Dynamic-Mesh-Repartitioning ..........................................................................73 6.2 Persistent Geometrical Induction matrix .....................................................................77 7. BOUNDARY PHYSICS ...................................................................................................79 7.1 Slip walls ......................................................................................................................79 7.2 Velocity inlets ..............................................................................................................80 7.3 Trailing edges...............................................................................................................83 7.4 Symmetry planes ..........................................................................................................87 8. RESULTS AND ANALYSIS ............................................................................................90 8.1 Basic geometry test suite .............................................................................................91 8.1.1 Rectangle wings ............................................................................................91 8.1.2 Delta wings ...................................................................................................93 8.1.3 Double-delta 75? / 65? ..................................................................................95 8.1.4 Arrow wing 70? ...........................................................................................97 8.1.5 Diamond wing 70? .......................................................................................99 8.1.6 Sweep wing 45? .........................................................................................101 8.1.7 Unit sphere and axisymmetric bodies .........................................................103 8.2 Advanced geometry test suite ....................................................................................106 8.2.1 F-16A ..........................................................................................................107 8.2.2 F-18A ..........................................................................................................114 8.2.3 DLR-F6 .......................................................................................................121 8.2.4 NASA Canard Fighter.................................................................................128 8.2.5 Rutan Varieze..............................................................................................136 9. CONCLUSIONS .............................................................................................................143 REFERENCES ......................................................................................................................148 APPENDIX .............................................................................................................................155 vii LIST OF ILLUSTRATIONS Figure 3.1: A typical structured mesh over a high curvature surface ......................................... 18 Figure 3.2: Solver run time data versus mesh size for various geometries................................. 19 Figure 3.3: Definition of manifold and non-manifold mesh surfaces ......................................... 21 Figure 3.4: Case study sphere in structured and unstructured environments ............................. 24 Figure 4.1.1: The mapped vortex ring for a facet on the surface mesh ...................................... 27 Figure 4.1.2: Velocity induced at a point by a segment of a facet-bound vortex ring ................ 28 Figure 4.3.1: The wake-strand model at the trailing edge .......................................................... 35 Figure 4.3.2: Vorticity decay rates and velocity desingularisation ............................................. 37 Figure 4.3.3: Adaptive pseudo-time stepping in the relaxed strand theory ................................ 39 Figure 4.3.4: Local physics driven wake strands ........................................................................ 40 Figure 4.3.5: Relaxed wake strands on bodies of revolution ...................................................... 41 Figure 4.4.1: Concept of integrated circulation .......................................................................... 42 Figure 4.4.2: Integrated circulation loops around a sample geometry ........................................ 43 viii Figure 4.4.3: The elliptical wing vorticity and wake solution .................................................... 45 Figure 4.4.4: Integrated circulation loops around the elliptical wing ......................................... 45 Figure 4.4.5: Composite span-wise sectional lift coefficient for the elliptical wing .................. 46 Figure 4.4.6: Integrated circulation loops around the F-18A geometry ..................................... 47 Figure 4.4.7: Composite span-wise sectional lift coefficient for the F-18A ............................... 48 Figure 4.4.8: Span-wise integrated circulation for rectangle wing (AR=4, 12? AOA) .............. 49 Figure 4.4.9: Lift coefficient as function of spacing between ICLs for rectangle wing ............. 49 Figure 4.4.10: Topological ?walking? used to determine B.L. growth length ........................... 51 Figure 4.4.11(a): Vorticity distribution for the Bell-Boeing tilt-rotor (5? AOA) ....................... 52 Figure 4.4.11(b): Skin-friction coefficient for the Bell-Boeing tilt-rotor (5? AOA) .................. 52 Figure 5.1.1: Size based meshing case study for a truncated cone surface ................................ 56 Figure 5.1.2(a): Initial pressure-solver mesh for the Varieze (108,040 facets) .......................... 58 Figure 5.1.2(b): Re-tessellated vorticity-Solver mesh for the Varieze (10,036 facets) .............. 58 Figure 5.1.3: Mesh sizing distribution and gradients .................................................................. 59 Figure 5.2.1: Unstructured mesh facet quality distributions ....................................................... 61 Figure 5.3.1: Edge collapsing on below-threshold facets. Before (left), after (right) ................. 62 Figure 5.3.2: Edge swapping on facet pairs, before (left) and after (right) ................................ 63 ix Figure 5.3.3: Individual facet-splitting, before (left) and after (right) ........................................ 64 Figure 5.3.4: Patch-splitting, before (left) and after (right) ........................................................ 65 Figure 5.3.5: Piercing correction, before and after (cutting surface hidden) .............................. 66 Figure 5.3.6: F-16A initial geometry pre-solver topological probing results ............................. 68 Figure 5.3.7: F-16A case study results for automated mesh refinements ................................... 68 Figure 5.4.1: Mig-21 mesh conformity case study results .......................................................... 71 Figure 5.4.2: Mig-21 case study convergence history results ..................................................... 72 Figure 6.1.1: Rationale behind FDMR ....................................................................................... 74 Figure 6.1.2: The F-35 geometry before and after FDMR ......................................................... 76 Figure 6.1.3: Convergence history for the F-35 solution (AOA = 5?, 1 m/sec flow) ................. 77 Figure 7.2.1: The F-18 geometry with velocity inlets and symmetry plane ............................... 82 Figure 7.2.2: The F-18 geometry with no Neumann relaxation ................................................. 82 Figure 7.2.3: The F-18 geometry with Neumann relaxation ...................................................... 82 Figure 7.3.1: Condition-1 for auto-detecting ?thick? trailing edges marked by the solver ........ 84 Figure 7.3.2: Centrifugal normal vectors for a general facet ...................................................... 85 Figure 7.3.3: Typical results of trailing edge auto-detection ...................................................... 85 Figure 7.3.4: The F-18 geometry manifold and non-manifold surfaces ..................................... 86 x Figure 7.4.1: The DLR-F6 geometry with and without symmetry BC ....................................... 88 Figure 7.4.2: The DLR-F6 geometry with symmetry BC (AOA = 0?, 1 m/sec flow) ................ 88 Figure 8.1.1: Lift coefficient versus angle of attack for rectangle wings ................................... 92 Figure 8.1.2: Lift versus Drag for rectangle wings ..................................................................... 93 Figure 8.1.3: Lift coefficient versus angle of attack for delta wings .......................................... 95 Figure 8.1.4: Unstructured mesh on the 75?/65? double-delta ................................................... 96 Figure 8.1.5: Lift coefficient versus angle-of-attack for 75?/65? double-delta .......................... 97 Figure 8.1.6: Unstructured mesh on the 70? arrow wing ............................................................ 98 Figure 8.1.7: Lift coefficient versus angle-of-attack for the 70? arrow wing ............................. 99 Figure 8.1.8: Unstructured mesh on the 70? diamond wing ....................................................... 98 Figure 8.1.9: Lift coefficient versus angle-of-attack for the 70? diamond wing ...................... 100 Figure 8.1.10: Lift coefficient versus angle-of-attack for the 45? sweep wing ........................ 102 Figure 8.1.11: Normalized sectional lift coefficient for the 45? sweep wing ........................... 103 Figure 8.1.12: The unit sphere .................................................................................................. 104 Figure 8.1.13: Lift coefficient versus angle of attack for the sphere ........................................ 104 Figure 8.1.14: The ogive axisymmetric body ........................................................................... 105 Figure 8.1.15: Lift coefficient versus angle of attack for the ogive body ................................. 106 xi Figure 8.2.1: The Lockheed-Martin F-16 ................................................................................. 107 Figure 8.2.2: Original geometry for the F-16A post CAD design ............................................ 108 Figure 8.2.3: Auto-mesher results on the F-16A ...................................................................... 110 Figure 8.2.4: Vorticity distribution on the F-16A ..................................................................... 111 Figure 8.2.5: ICL distribution about the symmetry plane on the F-16A .................................. 112 Figure 8.2.6: Lift coefficient span-wise distribution about the symmetry plane ...................... 112 Figure 8.2.7: Lift coefficient versus angle of attack for the F-16A .......................................... 113 Figure 8.2.8: Lift versus drag for the F-16A ............................................................................. 114 Figure 8.2.9: The McDonnell Douglas F-18 ............................................................................. 115 Figure 8.2.10: Auto-mesher results for the F-18A .................................................................... 115 Figure 8.2.11: Vorticity distribution on the F-18A ................................................................... 116 Figure 8.2.12: ICL distribution about the symmetry plane on the F-18A ................................ 117 Figure 8.2.13: Lift coefficient span-wise distribution about the symmetry plane .................... 118 Figure 8.2.14(a): Panel model of the McDonnell Douglas F-18A with 1,336 panels .............. 119 Figure 8.2.14(b): Unstructured model of the McDonnell Douglas F-18A with 7,853 facets ... 119 Figure 8.2.15: Lift coefficient versus angle of attack for the F-18A ........................................ 120 Figure 8.2.16: Lift versus drag for the F-18A ........................................................................... 121 xii Figure 8.2.17: The DLR-F6 geometry on a mount in the DLR wind-tunnel ............................ 122 Figure 8.2.18: Initial AIAA CAD tessellations for the DLR-F6 .............................................. 123 Figure 8.2.19: Auto-mesher results for the DLR-F6................................................................. 123 Figure 8.2.20: Vorticity distribution on the DLR-F6 ................................................................ 124 Figure 8.2.21: ICL distribution about the symmetry plane on the DLR-F6 ............................. 125 Figure 8.2.22: Lift coefficient span-wise distribution about the symmetry plane .................... 126 Figure 8.2.23: Lift coefficient versus angle of attack for the DLR-F6 ..................................... 127 Figure 8.2.24: Lift versus drag for the DLR-F6 ........................................................................ 128 Figure 8.2.25: The Canard Fighter model on a mount at NASA Ames tunnel ......................... 130 Figure 8.2.26: Auto-mesher results for the Canard Fighter post CAD design ......................... 130 Figure 8.2.27: Vorticity distribution on the Canard Fighter ..................................................... 132 Figure 8.2.28: ICL distribution about the symmetry plane ....................................................... 133 Figure 8.2.29: Lift coefficient span-wise distribution about the symmetry plane .................... 133 Figure 8.2.30: Lift coefficient versus angle of attack for the Canard Fighter .......................... 135 Figure 8.2.31: Lift versus angle of attack for the Canard Fighter ............................................. 136 Figure 8.2.32: The Rutan Varieze ............................................................................................. 137 Figure 8.2.33: Auto-mesher results for the Varieze post CAD design ..................................... 138 xiii Figure 8.2.34: Vorticity distribution on the Varieze ................................................................. 139 Figure 8.2.35: ICL distribution about the symmetry plane ....................................................... 140 Figure 8.2.36: Lift coefficient span-wise distribution about the symmetry plane .................... 141 Figure 8.2.37: Lift coefficient versus angle of attack for the Varieze ...................................... 142 xiv LIST OF TABLES Table 1.1: Overview of Aerodynamic Evaluation Schemes ......................................................... 2 Table 3.1: Case study results ...................................................................................................... 25 Table 5.1.1: Size based meshing results for the truncated cone surface ..................................... 56 Table 5.3.1: F-16A automated mesh refinement results ............................................................. 68 Table 7.4.1: DLR-F6 solver results comparison (AOA = 0?, 1 m/sec flow) .............................. 89 Table 8.2.1: Auto-mesher results for the F-16A ....................................................................... 110 xv NOMENCLATURE ? X-coordinate of spatial point P ? Y-coordinate of spatial point P ? Z-coordinate of spatial point P Velocity induced by the ? edge of a vortex ring mapped to a facet Circulation around a marked spatial loop ?? Incremental length along an edge of a mapped vortex ring Position vector of a spatial point P from a vortex filament Total edges on a polygonal facet of arbitrary valency Geometrical influence coefficient for a spatial point P by mesh facet ? ? Vorticity shed into the wake node ? Number of vorticity sources flowing into a wake node Number of vorticity sinks flowing away from a wake node Angle between the vortex filament length vector and the position vector from the terminal vertex of the filament to the spatial point P xvi Angle between the vortex filament length vector and the position vector from the initial vertex of the filament to the spatial point P (?) Radius of a vortex core along a given wake strand at pseudo-time ? ? Angular velocity of a two-dimension vortex ?? Increments in pseudo-time ? Distance to the downstream location of the Trefftz plane ? Distance to the seed point of the wake node ? Component of the free-stream velocity vector in a direction normal to the surface of the Trefftz plane ? Free-stream density ?? Incremental length along the integrated circulation loop span-wise direction ? Span of the integrated circulation loop surface around a geometry Free-stream velocity L Induced lift around a geometry ? Induced drag around a geometry ? Induced downwash between the two consecutive integrated circulation loops at the span-wise location ? ? Downwash evaluation position between any two consecutive integrated circulation loops xvii ? Free-stream fluid viscosity ?? Local Reynolds number for a given mesh facet ? Skin friction coefficient for a given mesh facet ? Length of edge between vertices 1 and 2 for a given mesh facet ? Length of edge between vertices 1 and 2 for a given mesh facet ? Length of edge between vertices 1 and 2 for a given mesh facet ? General vector potential field 1 1. INTRODUCTION Current industry trends point toward an increasing use of towards optimization based design pipelines for aerospace and marine applications 1-4. Numerous optimization schemes have been developed based on varying underlying philosophies 2-3. However, for complex systems, one aspect is common to all of these optimizers: the need for hundreds or thousands of design evaluations to arrive at a truly global optimum design 1-2. Since each of these designs must be tested for aerodynamic or fluid-dynamic performance, a balance between minimizing the computation time versus maintaining fidelity of solution and validity of the physics must be made. This balancing represents a major limitation in the comprehensive use of Multi-Disciplinary-Optimization (MDO) design pipelines 2. There is a certain subset of MDO applications in aerospace and marine designs where the determination of aerodynamic loads is the precursor for the quantified performance index in the optimization process 2. For these applications, traditionally four major methodologies have been used. These include empirical formulations, potential- flow methods, solution of Navier-Stokes or Euler equations in volumetric meshes (hereby generally referred to in this document as Computational-Fluid-Dynamics or CFD) and experiments in the field (wind- and water-tunnels as well as flight testing). Each of these four approaches brings with it a unique requirement for user inputs, computation times, 2 fidelity of solution, possibility for integration into MDO, range of applicable physics and overall cost of the process. Table-1.1 summarizes these parameters in a general sense for each of the four approaches. Clearly each approach balances out advantages in certain niche areas with disadvantages in other sectors. Table 1.1: Overview of Aerodynamic Evaluation Schemes Empirical equations Potential methods Navier-Stokes solvers Experimental methods Complexity of input data/geometry Low Intermediate Heavy Heavy Run time for single design Low Intermediate Intermediate/ Heavy Heavy Solution fidelity Intermediate Intermediate Good Good Automation for optimization Very high high Low Very low Range of physics Limited Limited High High Cost of process Very low Low High Very high Experimental methods usually offer the closest approximation to actual flight test data but at exceptional cost per design case. In addition, they are highly susceptible to environmental errors and anomalies (i.e. noise) in the results if any attempt is made to reduce the costs or the computational time for the test, all of which are the highest compared with the other three approaches. As such, experimental schemes are restricted 5 to studies of designs at advanced levels within the design cycle and typically follow CFD testing and validations of the design expectations. CFD approaches are rapidly making inroads into the MDO environment and have had visible success in demonstrating the use of wide range of physics during optimizations 2. However, the approach is currently restricted to simpler optimization cases for numerous reasons 2. One reason is the hardware associated with CFD approaches. For aerospace applications in the external-aerodynamics categories as well as certain marine applications, the necessary hardware to do thousands of CFD runs is very prohibitive to most users during the early design phases. However, this limitation is expected to be temporary. Current trends suggest that this limitation is rapidly being eroded as hardware expansion costs reduce industry-wide and shared industry pooling of large computing facilities becomes common practice 2. Even so, for the present time, it is expensive to use CFD in MDO pipelines. The primary limitation, however, is more philosophical in nature and relates to how CFD approaches are applied over a user specified geometrical surface. Since CFD approaches are volumetric in nature, i.e. they require volume meshes around the surface of interest to solve for the flow-field, the weak link in the solution fidelity chain becomes the quality of mesh and the automation of mesh generation. Current volumetric meshers are more robust today than at any other time as far as single run cases are concerned. However, they still require significant user inputs and visual verifications before being passed to the solver. For a large number of external aerodynamics problems, solver stability is directly related to the overall volume mesh quality while the solution fidelity 6 is dependent on the refinement of volume meshes in areas of interest around the input geometry. Slight variations in the mesh cell density or quality can cause solver divergence during iterations (non-physical ?trip-ups?) or cause fluctuations in the evaluated loads (from pressure field summations) that can render useless an automated MDO process. Furthermore, the requirement to visually identify regions of interest for further refinement prior to solver runs makes it virtually impossible to automate the process under the direction of an optimizer. Generally, the meshing is split into two phases in CFD approaches: surface and volume. The input geometry must be cleaned up in order to create a sufficiently refined surface mesh before the volume meshing can begin. This creates an additional layer of possible errors and user inputs and the final result is highly sensitive to the user settings during the surface meshing phase. Finally, there is the issue of interfacing with the optimizer. Because of the dozens of independent meshing parameters that must be varied for changing input geometries, automation of geometry variation is a major problem in terms of optimizer control 2-3. So far, only simple cases with modest design spaces have been shown to be optimized using CFD. This is expected to change in the coming years as more focus is applied on making the meshing independent of user inputs. Potential flow methods have their own unique set of advantages and limitations with respect to an MDO environment 5-7. Depending on the type of solver chosen for these approaches, the choice of user inputs, flexibility, stability and output data varies. 7 However, the biggest limitation of potential flow solvers is the restricted range of physics in which they operate 2, 5-7. This is typically the purely subsonic flow regimes (Mach < 0.8) or the truly supersonic flow regimes (Mach > 1.5) and even then only for inviscid flow-fields. By comparison, both CFD and experimental methods can smoothly transition between subsonic, supersonic and hypersonic flows with both viscous and non-viscous environments. However, for external aerodynamics MDO requirements for cruise conditions, this range of physics is usually very usable given the other major advantages offered by potential flow solvers. Typical applications in this genre include subsonic aircraft and missiles designed for long-range cruise flights, slow-speed wind-energy turbines, sail- boat designs, steady and unsteady flow behind propellers, aerostats and marine propulsion. Once the initial assumption is made to use potential-flow solvers inside their range of applicability, the other inherent advantages make them highly competitive and popular for MDO. One of the major advantages of potential-flow solvers over CFD solvers is the requirement of only surface meshes in order to capture the entire volumetric effect caused by the input geometry to the exposed flow-field. The elimination of volume meshing removes many of the aforementioned limitations regarding volume meshers. Furthermore, surface meshes result in substantial savings in computation time during the design evaluations (volume meshing for external aero cases can consume 20-40% of the overall solver time in CFD 3; this estimate applicable to most tetrahedral/polyhedral meshers, i.e. Delaunay meshers) as well as reducing the overall hardware requirements 8 (the majority of the memory requirements are associated with volume mesh generation and saving of the generated mesh on the local machines). Another advantage inherent to potential-flow solvers is the solver stability 2. Only geometrical surfaces with self-intersections, non-manifold interfaces and non-conformal patches on the surface mesh are known to be problem areas. Other areas of instability can include the inappropriate use of automated wake relaxations but this is generally minor. Finally, the solution convergence times associated with potential-flow solvers is known to be very efficient for most applications even when using surface meshes with large variations in mesh densities 2, 5. This reduces overall computation times for an MDO application substantially compared with a similar CFD approach. Potential-flow solvers are not without their inherent limitations beyond the restricted range of physics, however 5-7. As these solvers have become more advanced over the years, they have started to emulate the same limitations that afflict CFD solvers. The primary cause behind this is the use of pressure field summations to evaluate loads over the geometry. The need for this is driven by the attempts to push potential-flow solvers as alternates to CFD solvers for aerospace applications by various authors and organizations 5-11. For example, NASA?s PMARC, VSAERO 8-10 and PANAIR 11-12 potential flow solvers attempt to include schemes to evaluate viscous boundary layers 13 and flow separation 8-10 using the pressure fields evaluated by the underlying potential- flow solvers. This has been a semi-successful attempt to expand the physics of potential- flow solvers and move the methodology closer to CFD. However, doing so has meant 9 that these solvers have lost several of the aforementioned advantages associated with potential-flow solvers whilst also not attaining the true fidelity of the CFD solvers with regard to viscous flow regimes and flow separation effects. This argument can be elaborated on using a philosophical overview of the approaches taken by the three NASA potential-flow solvers. All three of these codes attempt to add a pseudo-coupled boundary layer model to the potential-flow results by adjusting the geometrical surface to account for the thickness of the local boundary layer 8-13. This boundary layer thickness is evaluated locally using the knowledge of the local pressure and velocity fields as well as the distance from the nearest leading edge marked by the pressure field results. And while they do gather some success 8-11 in establishing the boundary layer flow, the limitations incurred on the solver outweigh any substantive gains from such an approach. The primary limitation incurred is the requirement for these solvers to use a pressure-solver 8-11 rather than a vorticity-solver (elaborated later in this document). A pressure-solver necessitates the use of quadrilateral source distributions coupled with trailing edge doublet distributions 8-10. As such, the solver becomes highly sensitive to the fidelity of the surface mesh on which it is applied (similar to CFD solvers and unlike vorticity-solvers). The pressure-solver cannot tolerate local surface ?bumps? and ?dents? on the mesh without lowering fidelity of the evaluated pressure field and hence the aerodynamic loads 9-10. Furthermore the solver results become highly dependent on the accurate capture of the pressure field. This means that for simple airfoil cases, for example, highly refined surface meshes must be created in the high curvature regions in 10 order to capture the pressure gradients 11. This is again similar to the CFD solvers and causes similar cascading effects downstream into the optimization pipeline including higher computation times, solver sensitivity to surface meshes and reduced flexibility for the optimizer to dramatically vary the geometry within the design space 2. If further flow modifications are needed such as viscous boundary-layers 13, wherein the potential solution must be iterated until stability with the boundary layer is established, the computation times go up even higher 9-11, 13. Another attempt in recent times has been to incorporate separation flow effects within the potential-flow solutions 11. The general idea behind this is to use the pressure-field evaluated from the potential solver to determine flow separation points on the geometry using empirical formulations and use them to move the trailing edge locations up or down on the geometry. This is computationally very intensive and has been shown to be effective only under very controlled conditions 10-11. As such, any solution results for known flow separation cases are purely the domain of CFD solvers that can evaluate these regimes more efficiently and with relatively higher solution fidelity. Pressure-solvers and Vorticity-solvers use different panel types to suit their underlying theme. Pressure-solvers, as the name alludes, use the determination of the pressure field around given geometry to determine loads. As such, they use combination of source and sink panels coupled with doublet panels for the trailing edge. Source and sink distribution panels are mathematically accurate in their evaluation of the local pressures and velocities and hence suit pressure-solvers. 11 Vorticity-solvers take an entirely different approach altogether. They depend on the evaluation of the circulation around a given geometry and use the Kutta-Joukowski formulation applied to each panel for evaluation of the induced loads 14. As such, these solvers typically use vortex rings or doublet distributions over the surface facets 15. These types of potential flow panels are not suited for evaluation of the pressure field accurately except in far-field regions relative to the panel and on panel control points where the external physics conditions are satisfied 14. With these underlying philosophical differences, the requirement for surface meshes and their applicability to an MDO environment vary substantially between the two solvers 2, 14. Pressure-solvers are able to work with unstructured meshes as well as structured ones, mainly because each surface facet is evaluated for a certain source or sink strength regardless of the orientation of its edges relative to the flow direction. However, vorticity-solvers have thus far been used only with structured meshes since the local application of the Kutta-Joukowski formulation for loads at the facet level requires segregation of facet edges into two categories: bound or trailing 2, 14. This forces the mesh to rely on pre-established U-V mapping and restricts the use of advanced geometries with vorticity solvers 2. To balance this limitation, vorticity-solvers are known to be substantially more robust and stable compared to pressure-solvers 14. They allow for the user of non- manifold and non-conformal mesh surfaces (which the pressure-solvers cannot handle) and are less sensitive to surface perturbations. They also allow use of coarser meshes on 12 the geometry 14. This reduces the computation time substantially and affords greater stability to the solver, making it usable in an optimization environment 2. It is therefore apparent that both solvers have their inherent advantages and disadvantages. However, it is also apparent that in order to be truly applicable for a wide range of MDO applications, niche advantages from both solvers is needed and must be combined to harness the true potential of potential-flow solvers. This forms the basis of the current effort. 2. MOTIVATION A summary of the main objectives of the new approach can be written thus: A) Retain advantages of the potential-flow vorticity-solvers: ? Fidelity of solution: Aerodynamic loads (lift, induced drag and viscous drag) are used as a measure of the fidelity of the solution. The range of free-stream conditions is for Mach Numbers below 0.8 using the compressibility correction factor. It will be shown later that for certain geometries the range of free-stream Mach number can be extended to 0.9. ? Speed of solution: High solver speed is paramount for use of any aerodynamic prediction methodology within an MDO environment. It is possible to use the overall CPU time till solver convergence as a measure of the solver speed. This document will refer to solver speed in the context of this parameter throughout. ? Flexibility: Ability to evaluate aerodynamic loads over surfaces with high curvature, proximity, intersections, non-conformity of patch interfaces 14 and abrupt sizing gradients is paramount to creating a high fidelity MDO environment. Although difficult to quantify in and of itself, the measure of the speed of the solver and the fidelity of predicted loads will be used as a measure of the robustness against such ?high stress? geometrical surfaces ? Stability of solver: The solver must always attempt to prevent divergence during iterations even when applied over difficult-to-solve geometrical surfaces. Quantification of solver stability must be developed and used to improve automated prognosis methods to prevent solver divergence. B) Establish ability to work with unstructured surface meshes ? Allow use of triangular facets of arbitrary orientation relative to flow ? Bring potential-flow surface meshing to current industry standards: This can be split into several sub-goals: a) Introduce size based meshing to potential flow solvers b) Quantify surface mesh quality and automate its control c) Maintaining mesh intra-patch and inter-patch conformity or show lack of sensitivity of solver to it d) Introduce surface meshing data-packet structures similar to CAD/CAE industry standards 15 ? Remove requirement for forced U-V mappings on underlying geometry, thereby improving ? Open wide-spectrum geometries to potential-flow solvers. This remains a qualitative goal which can be demonstrated via a diverse testing suite for the current solver. Each of these goals will be addressed in the current document and a summary of the results will address the degree to which they have been achieved. 3. THE NEED FOR UNSTRUCTURED MESHES A structured surface mesh is typically defined as one where a given surface is mapped along two numerical axes (?i-j? or ?u-v?) and each facet is bound by four vertices extracted from the rectangular structure of the underlying map (known as the U- V map). In essence, this means that regardless of the physical convolution or warping of the physical surface, the underlying U-V mapping is rectangular with uniformly spaced increments. Each vertex of this mesh has an associated pair of numerical coordinates so that a vertex P is stored in the mesh as: (3.1) And a surface facet is created out of a set of four corner vertices as: (3.2) As such, any structured surface mesh facet is essentially quadrilateral in nature. As a result, facets have to be ?forced? to assume triangular shapes to meet geometrical demands. This is accomplished by ?edge-collapsing? wherein two of the facets on a given ?? = [?(? ?) ?(? +1 ?) ?(? +1 ? +1) ?(? ? +1)] ?? = [?(? ?) ?(? +1 ?) Y(? +1 ? + 1) ?(? ? +1)] ?? = [?(? ?) ?(? +1 ?) Z(? +1 ? +1) ?(? ? +1)] ? = ?(? ?) ? = ?(? ?) ? = ?(? ?) 17 edge of the facet are merged on top of each other to reduce the quadrilateral facet to a triangle in physical space. In numerical space such triangles still retain the underlying quadrilateral U-V map. The advantages and disadvantages of such approaches to mesh generation are obvious. The forced quadrilateral U-V mappings of the underlying facets ensure that the memory requirement for a given forced structured triangle is always more than an equivalent unstructured triangle facet. In addition, forced U-V mappings ensure that regions of high curvature or concave warping are always populated with more facets than actually required from a geometrical standpoint. As a result, a structured mesh patch forced on such geometry always results in higher mesh facets. Further, the quality of such facets is typically very poor (with downstream in the solver where oscillations are induced). Additionally, the typical quadrilateral facet has to be sanitized to establish its ?effective? surface normal direction given its higher level nature with regard to the baseline triangle. With triangular facets, this is not required given the deterministic nature of the triangle geometry in three-dimensional space. However, since the quadrilateral facet is a combination of a pair of triangles (which can be further varied depending on choice of diagonal used), there is no fixed surface normal. This problem is typically eliminated by evaluating the surface normal as a cross-product unit vector result of the two diagonal vectors of the facet. Understandably, this reduces the geometry capture fidelity of the facet. If the quadrilateral facet were to be split into two facets to improve this fidelity, the number of resultant mesh facets increases. 18 Figure-3.1 illustrates the application of a structured mesh patch on a high curvature surface. The long needle-shaped facets near the nose of the geometry are visible. Such facets are typically of very poor quality and can significantly increase the mesh facet count. Figure 3.1: A typical structured mesh over a high curvature surface A general rule of thumb obtained from the current effort shows that the computation time for the solution convergence is directly proportional to the square of the mesh facet count (Figure-3.2). ????????????? ???? ? ????? (3.3) 19 Figure 3.2: Solver run time data versus mesh size for various geometries. All geometries run for 10 iterations only Further, the regions of lower curvature are populated with much higher number of facets than needed because of the need to satisfy the underlying U-V mapping with appropriate physical locations. Such adherence to the U-V mapping means that the mesh is unable to recognize and react to high curvature regions effectively in physical space. This restricts the effective application of these meshes to simplified regions with long gradual curvature regions. Despite the aforementioned shortcomings of the structured surface meshes, with potential-flow solvers (and especially the vorticity-solvers) they offer certain unique advantages. Arguably, vorticity-solvers in their prior form have been made possible mainly because of the inherent advantages of the structured mesh. Further, mesh 20 generation with structured surface mesher is much simpler from a development standpoint than a comparative unstructured surface mesher. The vorticity-solvers require a facet edge pool to be split into bound and trailing vortices during the application of the Kutta-Joukowski formulation for induced loads at the facet level. This means that facets of arbitrary orientation (as obtained from unstructured meshes) cannot be applied. Traditionally, the U-V mapping of the structured meshes has been used in all legacy vorticity-solvers to establish the bound and trailing vortices. The bound vortices are simply marked as those facet edges aligned along a certain numerical axis (U or V) and the other edges are marked as trailing. Note that such marking usually only corresponds to the correct flow field if the physical space warping and curvature is not significantly different from the numerical axes alignment. Therefore, such solvers are only applicable on geometrical surfaces that meet this criterion which is a major limitation of the vorticity-solvers. Further, the induced load formulation requires the knowledge of the vorticity in the next bound vortex downstream of a given facet. This information is also evaluated from the U-V mapping information. This reduces the computational burden on the solver and allows the evaluation of loads in a simple and straightforward manner. Having established the merits and demerits of the structured surface meshes, a similar analysis can be done for the unstructured surface meshes. The primary philosophical difference between the two meshes is the lack of an underlying U-V map on the unstructured meshes and that the unstructured meshes are first and foremost 21 triangulated (instead of quadrilateral in structured meshes). Meshes of this type are not bound by any numerical forces in the distribution of surface vertices. Each surface mesh generated using this approach is essentially a numerical pool of vertices, edges and faces that are marked using connectivity mappings evaluated at the time of mesh generation. A given vertex may be shared by numerous edges depending on the curvature of the region. A given edge may be shared by two faces (for a manifold surface) or more than two faces (for a non-manifold surface) as shown in Figure 3.3. Figure 3.3: Definition of manifold and non-manifold mesh surfaces (selected facets are shown in pink; note the number of faces sharing an edge) Information storage for the mesh can vary. Industry practitioners use differing mesh storage concepts. For example, a facet can be fully defined by knowing its three vertices in three-dimensional space. Such a mesh representation can be written as F?V where F=Facets and V=Vertices. Similarly, a facet can also be defined entirely from the knowledge of its three edges (which in turn are marked from the vertex mapping) and such a mesh representation can be written as F?E?V . For the current solver the mesh information storage is F?V . 22 Data packets store each facet?s information including vertices, edges, topological connectivity and external physics as primary data while derivative data such as facet area, quality metrics, surface normal vector, centroid etc. are stored or evaluated as necessary. Unstructured surface meshes bypass some of the aforementioned limitations of the structured surface meshes. Geometry capture with unstructured meshes is more efficient and can be designed to significantly reduce overall mesh count and improve the quality of the facets (i.e. reduce the presence of large number of needle shaped facets). They can capture regions of concave and convex warping with much more even distribution of facets and facet edge-sizes. They can also be used to enact local surface refinements more effectively in physical space rather than having to do so in the underlying U-V space. Because each facet on such a mesh is a self-contained entity, bad quality facets can be deleted, improved or refined to meet mesh requirements. Refinement can reduce the overall solver convergence oscillations (caused by the poor-quality cells; discussed later) and thereby reduce the computation time to convergence. Unstructured meshes also allow superior control on the mesh facet sizing as well as the sizing gradients. Further, each facet can be marked with its local external physics and boundary information. This allows much simplified creation of boundary zones (such as velocity inlets and trailing edges), but also assists in the automated detection of trailing edges on a surface. The limitations of the unstructured surface meshes pertain mainly to their application with the potential-flow vorticity-solvers. And while the local evaluation of the 23 vorticity can be made similar to the pressure-solvers using the external and boundary physics on the surface mesh, the induced loads cannot be evaluated without marking appropriate bound and trailing vortices as discussed for structured meshes. This has been the main limitation in the use of unstructured meshes with vorticity solvers thus far and is the primary motivation for the development of Method of Integrated circulation as will be discussed in the following Sections of this document. However, before moving to the discussion of the solver and evaluation of integrated loads, a case study has been presented to further illustrate the need for unstructured meshes. A simple geometry has been taken (a unit radius sphere with its center at the coordinate axes origin) and meshed using a structured and unstructured surface mesher developed for this effort. The two spheres have been refined to similar mesh requirements, i.e. the number of mesh vertices along any quadrant arc of the sphere must be equal to ten (this can be visually inspected from Figure 3.4). Note that the structured mesh is forced to convert its quadrilateral facets into triangular ones so that they can both be run through the current solver. Additionally, note that the forced U-V mapping on the structured mesh causes the existence of a large number of bad quality facets near the polar regions of the sphere (Figure-3.4) whereas the unstructured mesh easily bypasses this potential pitfall. As a result, the final mesh size between the two meshes is substantially different (Table 3.1). The two meshes were run through the solver under identical conditions and the results are shown in Table-3.1. Note that the aim here was to obtain a potential flow solution around the two spheres. As a result, no wake definitions were necessary and the 24 expected coefficient of lift for any angle-of-attack was expected to be zero. The flow velocity for the current case study was kept at 1.0 m/sec. Figure 3.4: Case study sphere in structured and unstructured environments The unstructured surface mesh is seen to take about ~88% less time to reach much superior convergence limits in fewer iterations and with higher solution fidelity. The structured mesh reaches the limit set on the solver iterations and still does not meet the required convergence criteria. This is caused by the solver oscillations resulting from large number of poor quality facets in its polar regions. The resulting mesh size created out of a similar refinement criterion is also favorable for the unstructured mesh (648 facets) compared with the structured mesh (1,404 facets). The resulting trends of this case study are applicable to a wide range of test cases and geometrical variations and help establish the necessity and advantage of the unstructured meshes for the vorticity-solvers. The next Sections of this document discuss the development of the Method of Integrated Circulation to allow evaluation of circulation based loads on such meshes. 25 Table 3.1 Case study results Structured Unstructured Reduction Number of facets 1404 648 53.85 % Solver run time (Seconds) 148.88 17.59 88.19 % CL (expected = 0.0) 0.000025 0.000009 64.00 % Convergence (expected = 1E-3) 9E-3 7E-4 92.22 % Solver iterations (max = 20) 20 7 65.00 % 26 4. SOLVER METHODOLOGY This Section details the vorticity solver used for the current effort along with the modified approaches for the evaluation of the integrated aerodynamic loads. The Sections are split to include discussions on the following: a) The generalized spatial velocity inductions by segmented elements of the vortex ring mapped to each facet b) The philosophy of the solver from a numerical application standpoint c) Development of the relaxed wake-strand model for unstructured meshes d) The method of integrated circulation for evaluation of the induced loads on the unstructured mesh as well as the modifications to skin-friction drag models to allow their application into a vorticity-solver 4.1 SPATIAL VELOCITY INDUCTIONS At its most fundamental level, the solver developed for the current effort reduces to the application of a vortex ring on a triangulated facet. Each edge of the facet therefore becomes a segmented part of the continuous vortex ring on that facet. Since each of these edges remains a finite, line segment in three-dimensional space, the induced velocity of 27 each segment at any other point in space is easily determined from the geometry of the plane described by the two end vertices of the line and the point of induction in space. Once this induction has been evaluated for each edge of the facet, they are summed up and treated as the induced velocity of the facet at the point of interest in space. Additionally, the direction of vortex ring is aligned to match the winding of the facet as shown in Figure-4.1.1. Because each valid (non-degenerate) facet on the surface mesh is a closed loop with a known winding, the vortex rings are easily mapped on their surface. This is also in accordance with Helmholtz laws which require every vortex ring to be closed unless stretched to infinity on both ends of a segment 14-15. Figure 4.1.1: The mapped vortex ring for a facet on the surface mesh The derivation for the velocity induced by a linear segment of a facet can be mathematically derived in terms of the coordinates of the two end vertices and the point 28 of interest. Consider a point P in space where a segment of a vortex ring is inducing a velocity as shown in Figure-4.1.2. Figure 4.1.2: Velocity induced at a point by a segment of a facet-bound vortex ring The velocity induced at P by an infinitely small segment of the ring edge is known from the application of the Biot-Savart law 14 as: ? = 4 ?? (4.1.1) Considering the geometry shown in Figure-4.2, Equation-4.1.1 can be modified: ? = 4 ?? ??( ) (4.1.2) Now, given that, ? = ?? ( ) 29 ? = ????( ) ?? = ? ? ( ( )) (4.1.3) Substituting Equation-4.1.3 into Equation 4.1.2, we get: ? = 4 ( ( )) ? ( ( )) ??( )? = 4 ??( )? ? (4.1.6) Equation-4.1.6 can be integrated along the entire length of the segment so that the formulation for the net induced velocity by one segment of the vortex ring is: = ? ? = 4 ?? ??( )? = 4 ?( ( ) ( )) (4.1.7) In Equation-4.1.7, all variables on the RHS are easily evaluated from vector formulations of the planar geometry defined by the three points in space (two vertices and one point of interest). Further, if the point of interest were to be so located that the radius becomes zero (all three points are collinear) then the solver inputs a zero velocity induction even though from Equation-4.1.7 it is obvious that we would have a mathematical singularity. Note that Equation-4.1.7 also reduces to the two-dimensional vortex Equation for the limiting case of infinitely long vortex (as a result of this, the two angles in Equation-4.1.7 becomes 0? and 180? respectively). Note also that the velocity magnitude evaluated by equation-4.1.7 is in a plane perpendicular to that defined by the 30 three points in space. However, knowing the surface normal vector for this plane from the three points, the velocity magnitude is broken down into three vector components along the global mesh coordinate axes. Equation-4.1.7 is the induced velocity by one segment of a facet. Therefore, to evaluate the velocity induced by the entire ring, we evaluate Equation-4.1.8 over all the relevant edges of the ring to get the net velocity induction as: = ? (4.1.8) A discussion of the derivation for the Biot-Savart law in the form used in Equation-4.1.1 is given in appendix-A. 4.2 SOLVER PHILOSOPHY Considering Equation-4.1.7, it is clear that the only unknown variable in the Equation is the vorticity strength. The LHS of the Equation is a variable that can be controlled at certain points on the geometry using the boundary physics conditions. The remaining terms in the RHS minus the vorticity strength are geometrically dependent and therefore determined from processing the surface mesh. If the vorticity strength could be determined using known values of the velocity at the control points, then the flow-field can be determined at all points in space around the geometry, allowing for evaluation of 31 the loads. This can be further quantified by considering Equation-4.1.7 and Equation- 4.1.8 as: = ? [ 14 ? ( ( ) ( ))] = (4.2.1) The coefficient is the geometrical influence coefficient for the velocity induction at point P by the facet ?j? of the surface mesh. If the point P were a control point, then the LHS is a boundary defined velocity term. Conceptually expanding Equation-4.2.1 for a mesh with N facets, the overview of the fundamental solver philosophy is: (4.2.2) Where the matrix B is the velocity requirements on the control points (and hence known from the external physics requirements) and the matrix A is the N?N matrix of geometrical influence coefficients for each of the control points. The matrix ? is therefore the N unknown variables that must be solved for the flow-field to be fully defined around the geometry. The system defined in Equation-4.2.2 is essentially a system of N Equations and N unknowns that can be solved explicitly. [ ? ? ? ? ? ? ? ? ? ][ ? ] = [ ? ? ? ] 32 The control points selected for triangular facets lie on the facet area-centroid and are unique in space (they have no edges from any facet on the mesh that crossed through it to introduce singularity) for a manifold, conformal surface mesh. The N?N matrix formulation in Equation-4.2.2 introduces certain problems as the size of the mesh increases. Since these matrices must be inverted during the solver run, the issue of numerical errors cascading up the inversion loop is severe and can lead to completely non-physical and erroneous results for the vorticity 2. Consequently, the matrix must either be evaluated in patches numerically during the inversion phase or the solver must be handed small patches of the mesh individually 16. The current effort handles this at the physical level and splits the mesh into partitions that are fed to the solver individually, freezing the remaining Sections of the mesh (their vorticity values) at the previous iteration. This approach creates an iterative solver for the overall surface vorticity. Given the nature of the geometry, an exact solution can never be arrived at practically and the solver iterations must be terminated once the surface vorticity values have ?converged? to a reasonable level of numerical accuracy. The convergence criterion used in the current effort is based on the stability of the surface vorticity solution and the stability of the relaxed wake strands in the Trefftz-plane downstream of the geometry. The vorticity convergence is evaluated as an area-averaged value of the facet vorticity strengths on the entire surface mesh. The wake-instability is quantified as the 33 RMS change of the two-dimensional strand positions within the Trefftz plane (Equation- 4.2.3). ????? ????? = ? [(? + ? ) (? + ? )] (4.2.3) The wake-instability criterion is secondary in nature and is under user toggle control. This is because terminating wake relaxation at a lower convergence value can lead to a wake that is not force-free and therefore can adversely affect the solution. Only in the case of simple geometries can such relaxation be frozen without affecting the overall result. Wake instability as a convergence criterion comes into play only as a means to reduce overall computation time and is entirely case dependent. A further note regarding compressibility can be made at this point. The application of the Biot-Savart law for velocity induction, although formulated for incompressible flows, can be modified with the addition of the compressibility correction factor 24 and is used in the current approach to simulate higher Mach number flows. The importance of such corrections is made clear in Section-8 of this document. 4.3 VORTICITY WAKES: THE RELAXED STRAND MODEL Unstructured meshes drive the wake geometry to be unstructured since the trailing edges are not mapped numerically. As a result, the wake sheet concept 2, 14 used in the 34 standard relaxed wake models is not applicable with the current approach. This problem has been resolved by the evolutionary extension of the relaxed wake sheet into the relaxed wake strand approach developed for this effort. The general principle for this approach is based around the conservation of vorticity 14 in the wake (in accordance with Helmholtz Laws of vorticity) for a strand emanating from ?nodes? rather than edges on the geometry. This is illustrated in Figure- 4.3.1. The starting point for this model is based on user-marked trailing edges on the geometry (either automatically or manually; see Section-7 on boundary conditions and external physics). Once the trailing edges are marked, the solver determines the vertex pool corresponding to these edges from the global vertex map. These vertices become the starting point or nodes, for the individual wake strands. Each node is then numerically propagated downstream using the induction effects of the overall facet-bound vortex rings, the free-stream velocity and the wake-strands from the previous iteration (excluding the strand being updated). The evaluation of the strength of the vorticity strand being propagated into the wake requires the knowledge of contributing ?source? and ?sink? edges for each node. These are obtained from the global vertex and edge connectivity maps and are illustrated visually in Figure-4.3.1. Edges of any facet that share the wake node are automatically marked as contributors to the node (unless the edge is marked as a trailing edge, in which case no vorticity is allocated to it and hence cannot be a contributor; green-highlighted edges in Figure-4.3.1). A vorticity contributor can be a source or a sink depending on 35 whether the direction of ?flow? of vorticity along the edge (determined from the facet winding) is towards or away from the wake node. Consider the two cases shown in Figure-4.3.1. Here the case on the left has three edges that share the wake node in addition to the two trailing edges also sharing it (marked green). The three edges sharing the wake are in turn shared by two facets each as part of the manifold surface. As a result, each edge is contributing the difference between the vorticity of its shared facets to the wake node. So this wake node has three contributing edges, two trailing edges and four contributing facets on the top surface. A similar situation may exist on the bottom surface as well. Figure 4.3.1: The wake-strand model at the trailing edge Similarly, consider the second case on the right in Figure-4.3.1. Here the wake node is being shared by two facets on the top surface, has two contributing edges and one trailing edge. The net vorticity of this wake then is the difference between the facet vorticity of the two facets sharing it. 36 Mathematically, the vorticity of the strand can be written as: (4.3.1) Because no new sources or sinks are added by the wake nodes, the net vorticity of the geometry remains conserved and the wake propagates only the required vorticity in the wake to satisfy the Kutta condition at the trailing edges. Vorticity strength reductions downstream of the trailing edge can only occur due to viscous effects. This is modeled using the application of the Lamb-Oseen 17-19 vortex model obtained from the exact solution of the Navier-Stokes equations for a laminar two- dimensional vortex. This model assumes a two-dimensional vortex with circular symmetry in which the streamlines are circles around the axis and the vorticity is a function of the radial distance away from the filament axis and the pseudo-time 17. The exact solution of the Navier-stokes equations 17 then takes the form: ? = ?? 4 ???( ) (4.3.2) Assuming the knowledge of the initial vorticity shed into the strand at the node on the trailing edge, Equation-4.3.2 is modified to give us the vorticity decay equation 17: ? = ? (1 ?( )) (4.3.3) ? = ? ? + ? ? 37 In Equation-4.3.3, the time used corresponds to the pseudo-time and starts at zero at the trailing edge. It is marched alongside the strand until the termination of the wake downstream at the Trefftz location. The behavior of Equation-4.3.3 is shown in Figure- 4.3.2 for the vorticity decay rate versus wake pseudo-time for a fixed radial position and for the velocity induction by a single filament in the radial direction for a fixed pseudo- time. The advantage posed by the implementation of the Lamb-Oseen model is visible from the results in Figure-4.3.2. The vorticity decay rate results in the ?desingularisation? of the rectilinear line vortex, in which the vorticity has a delta-function singularity. The velocity function also indicates the growth and an effective solid-rotation core along the length of the propagating wake-strand filament which simulates viscous behavior. Figure 4.3.2: Vorticity decay rates and velocity desingularisation using the Lamb- Oseen vortex core model 38 The numerical propagation of the wake nodes from the trailing edges downstream is evaluated using this pseudo-time increment evaluated for each step of each strand. The unstructured nature of the strand-based wake allows for the evaluation of strand-specific pseudo-time marching. This allows for the local ?compression? and ?stretching? of individual strands so that all of them are propagated to the same Trefftz plane location downstream. The Trefftz plane location remains the only external user specification for the generation of the relaxed wake strands. Once the time Trefftz location is known and the starting node on the trailing edge is evaluated, the local steps from any node location at a given time is evaluated using Equation-4.3.4: (4.3.4) This allows for the time step for each propagation step to change depending on where the previous propagation has taken the node relative to the fixed Trefftz plane location and the local velocity conditions. The final propagation step is either truncated or projected to intersect with the Trefftz plane. The time increment is thus adaptive and responds to each strand node individually. This is an advantage over wake-sheet methods since strands starting from differing positions on the geometry can be made to propagate to the same location downstream. This is illustrated in Figure-4.3.3. ?? = *? ? ? + ? 39 Figure 4.3.3: Adaptive pseudo-time stepping in the relaxed strand theory The wake strand model has several inherent advantages. The presence of adaptive time stepping allows for the elimination of all user inputs other than the downstream location of the Trefftz plane. User-specified time steps are no longer required. Further, the model allows the wake to conform to the geometry based on the local conditions. This is illustrated in Figure-4.3.4 for the case of the Rutan Varieze. The presence of the canard for this geometry creates relatively complex wake physics where certain strands follow their local physics and go above the wing whereas in other cases the go under the wing. Because each strand is a self-contained vortical entity, the wake-strand model allows such behavior to take place without problems. 40 Figure 4.3.4: Local physics driven wake strands Another advantage allows for by the strand model is the ability of the wake to be applied to bodies of revolution. This ability derives from the nodes based approach to vorticity shedding instead of the edge based models of prior wake-sheet models. A body of revolution has a trailing edge that rolls up on itself with the final node touching coincident in space with the first node. Since the wake-strand model makes no distinction as to the identity of vorticity contributing elements, the geometric relevance is unimportant. As a result, wake stability remains unaffected. This is particularly useful for geometries with podded engines as shown in Figure- 4.3.5 for the DLR-F6 geometry. Note the ?pull? generated by the wing-bound vorticity on the engine wake as it propagates downstream. It is possible that future enhancements can be made in conjunction with volumetric sources to allow study of integrated aero- propulsive effects on wing design. 41 Figure 4.3.5: Relaxed wake strands on bodies of revolution 4.4 AERODYNAMIC LOADS Aerodynamic loads are evaluated as three components in the current approach: a) The induced lift from the Method of Integrated Circulation b) The induced drag from the Method of Integrated Circulation c) The skin-friction drag from the local vorticity distribution 4.4.1 METHOD OF INTEGRATED CIRCULATION At the core of this approach is the extension of the basic laws on which potential- flow methods are based: linearization of aerodynamic effects 14. The extension is easily 42 illustrated in two dimensions as shown in Figure-4.4.1. The net circulation in the large loop encompassing the series of smaller loops with their own circulation is evaluated as: ? = ? + ? + ? ? + ? (4.4.1.1) Figure 4.4.1: Concept of integrated circulation This idea can be extended to three dimensions for a body of arbitrary surface tessellation. Assume that such a body has been solved for the local facet-bound vorticity using the vorticity solver described previously. Given that the tessellations are arbitrary, no application of the Kutta-Joukowski lift and induced drag equation can be made at the facets on this body, even though the vorticity values are evaluated correctly. However, given the potential-flow environment around this body, the net spatial effect of all the facets on the body at any point in space is equal to the integrated effect of the individual facets in isolation. If a two-dimensional closed loop were created around this body in a plane perpendicular to the direction of flow, we could evaluate the net 43 circulation around this loop by determining the velocity at all of the loop?s edges and then use the generalized equation for circulation 14 as: ? = ? ?? (4.4.1.2) If this were a two-dimensional body, the net induced lift from the body could be evaluated using the Kutta-Joukowski theorem 14. However, the body is three dimensional. As a result, several more sections must be created until the body is enclosed inside an effective bounding box of sectional planes, all of which are perpendicular to the direction of flow. An example of this approach is shown in Figure-4.4.2 for the Bell-Boeing advanced tilt-rotor geometry. Each of these loops is referred to in this document as an Integrated Circulation Loop (ICL). Once these loops have been created around the geometry of interest, the Kutta-Joukowski equation for induced lift is applied for each loop and the lift is evaluated for each slice of the bounding box in the span-wise direction. Figure 4.4.2: Integrated circulation loops around a sample geometry 44 The presence of these closely spaced ICLs also allows for the evaluation of the composite span-wise circulation distribution and is then used to determine the induced downwash between the ICLs using standard lifting-line equations 14 and using the net circulation value of each loop as the ?bound? circulation values. The induced lift and drag are thus evaluated as: L = ? ? ? ?? (4.4.1.3) ? = ? ? ? ? ?? (4.4.1.4) ? = 14 ? ( 1? ? ) ?? ?? ?? (4.4.1.5) Examples of the composite span-wise distribution function are given for two geometries: an elliptical wing and the Boeing F-18A Hornet in order to illustrate the nature of the function and its relation to the underlying geometry. Figure-4.4.3 illustrates the results for the elliptical wing of aspect ratio 5.1 tested at 10? AOA. The composite span-wise sectional-lift distribution is shown in Figure-4.4.5(a). Given the elliptical nature of the wing plan-form, the expected nature of the circulation distribution is elliptical 20 in nature and matches the predicted data within numerical accuracy limits (Figure-4.4.5(b)). The current approach reduces to the concept of the lifting-line solution for flat wings in the limiting condition. 45 Figure 4.4.3: The elliptical wing vorticity and wake solution Figure 4.4.4: Integrated circulation loops around the elliptical wing 46 Figure 4.4.5(a): Composite span-wise sectional lift coefficient for the elliptical wing Figure 4.4.5(b): Composite span-wise sectional lift coefficient for the elliptical wing 47 Figure-4.4.6 illustrates the same process applied for the Boeing F-18A geometry. The geometry was evaluated at 1? AOA with a symmetry plane and velocity inlets for the engine intakes under the wing. The auto-mesher was used resulting in a fine mesh and smaller average edge-size. This led to the solver using 44 ICLs in semi-span direction as shown in Figure-4.4.6. The complexity of the combined main wing, stabilizers and LERXs is captured in the span-wise lift distribution shown in Figure-4.4.7. Note that the spikes and valleys in the distribution function correspond to the additive or subtractive effects of the horizontal and vertical stabilizers, the LERXs and then the fuselage near the span-wise zero location in Figure-4.4.7. Figure 4.4.6: Integrated circulation loops around the F-18A geometry 48 Figure 4.4.7: Composite span-wise sectional lift coefficient for the F-18A In numerical application, the ICLs are discretized using the average edge-size of the underlying surface mesh. The spacing between ICLs is also discretized using the same reference values. It is possible to include other discretization values and this can have a varying effect on the solution fidelity depending on the complexity of the geometry in question. An example of a case-study is presented for a rectangle wing (AR=4) for 12? AOA. The spacing was varied about the average edge-size by different multipliers and the overall coefficient of lift evaluated for each of these distributions. The results are plotted in Figure-4.4.8 and Figure-4.4.9. It can be seen from the figure that increasing the spacing between the ICLs beyond the average edge-size gradually reduces the fidelity of the solution where reducing the distance between ICLs (and hence increasing the total number of ICLs) improves the solution fidelity slightly. 49 Figure 4.4.8: Span-wise integrated circulation for rectangle wing (AR=4, 12? AOA) Figure 4.4.9: Lift coefficient as function of increasing spacing between ICLs for the rectangle wing (AR=4, 12? AOA) 50 4.4.2 SKIN-FRICTION DRAG EVALUATIONS Evaluation of parasite drag in legacy potential-flow solvers has been based on the evaluation of the local Reynolds number for each facet 2, 21-22. This local Reynolds number is then used as part of the semi-empirical skin-friction coefficient equations. The most popular of these equations for high speed subsonic aircraft was developed by Prandtl-Schlichting 21-22 in 1932 and was subsequently modified with the turbulent correction factor. This equation can be written as: ? = 0 455 (l g ?? ) 1700?? (4.4.2.1) The application of this equation in a pressure-solver was made possible by the combined use of structured meshes and pressure-field evaluations for calculation of aerodynamic loads 8-11. The knowledge of local velocity and distance from the leading edge allowed local evaluation of the skin-friction coefficient for each facet that was then summed up for the net integrated value 8-11. With the current approach, the local velocity is not determined. However, the local distance from the leading edge can be determined by topological ?walking? for each facet towards the closest leading edge in the direction of the free-stream (Figure-4.4.10). 51 Figure 4.4.10: Topological ?walking? used to determine B.L. growth lengths Additionally, the local vorticity is known in lieu of the local velocity. As such, the equation for the local Reynolds number is heuristically modified thus: ?? = ? (? + ? ) ? (4.4.2.2) This model is applied after the determination of the local facet vorticity following convergence of the solution. The model follows certain heuristic arguments that were used in its formulation: a) When the facet is parallel to the flow, the vorticity on it is zero and equation 4.4.2.2 reduces to the standard formulation for a flat plate. 52 b) When the facet is at the leading edge, the local vorticity dominates the skin- friction evaluation since the distance from the leading edge is a very small value c) The dependency of the local Reynolds number on the surface vorticity means that the distribution patterns for the skin-friction coefficient closely resemble the vorticity distributions. This is also expected from theoretical arguments. An example of this is shown in Figure-4.4.11 for the case of the Bell-Boeing Advanced Tilt-rotor geometry. Figure 4.4.11(a): Vorticity distribution for the Bell-Boeing tilt-rotor (5? AOA) Figure 4.4.11(b): Skin-friction coefficient for the Bell-Boeing tilt-rotor (5? AOA) 53 5. MESHER METHODOLOGY This Section outlines the automation adopted for the unstructured surface mesher developed for this effort. The quantification of various heuristics and empirical parameters that control the overall surface mesh and help convert an input geometry into a solver ready surface mesh are also discussed. Finally, the Section will outline the importance of sizing and quality control on the solver mesh and the inherent robustness of the current approach to stressful input geometries. 5.1 SIZE BASED MESHING One of the opportunities provided by unstructured meshes is the ability to control sizing and sizing gradients for the surface facets at the local and global levels. This is of great importance in CFD simulations because of the volumetric cell based solver physics. In potential-flow solvers, the importance of mesh facet sizing varies depending on the choice of solver. Pressure-solvers are significantly more sensitive to the sizing and sizing gradients very similar to CFD solvers. This is not surprising in lieu of the similar nature of their aerodynamic load evaluators. As mentioned in previous Sections, capturing the pressure field accurately is of paramount importance for such solvers and this is not 54 possible unless the mesh is appropriately refined based on the local surface curvatures. Furthermore, surface patches making up the geometry must be conformal along the interface edges. This is due to the greater sensitivity of such solvers to gaps or punctures in the surface topology. If the solver is fitted with a boundary layer model, no such gaps can be allowed and the mesh surface across patches must be conformal at all times. Such conformity and sizing can be enforced through the careful development of an underlying U-V mapped structured mesh and is the reason for widespread use of such surface meshes with legacy potential-flow pressure-solvers. Vorticity-solvers are much more resilient to such limitations and afford greater flexibility in the sizing and sizing gradients. Effective solutions can be obtained even from meshes with sharply varying sizing gradients along the topology. Mesh conformity is not a necessary condition for solution fidelity. However, maintaining conformity along patch interfaces and a uniform sizing gradient greatly helps in faster solution convergence (boosted by lower oscillations in surface vorticity along solver partition interfaces as also because of higher quality facets). The savings in computation time due to sized based meshing therefore translate directly to faster runs within an MDO environment and therefore merit further investigation. Despite these advantages, there are certain pitfalls that must also be investigated. The primary disadvantage of size based meshing for vorticity-solvers is that it can sometimes defeat the very advantage that makes this approach MDO friendly. This can happen by the liberal use of size-based meshing over the entire surface geometry, which can significantly increase facet count in regions of low curvature where the 55 vorticity gradients are low and insignificant. Surface re-meshing, as brought about by the use of size based meshing is also particularly invasive to the smoothness of the initial surfaces and can create ?bumps?. However the vorticity-solvers are generally insensitive to local surface bumps while the Method of Integrated Circulation concept allows it to smooth out the effects of such topological perturbations. Finally, unless care is taken and appropriate metrics used, redefining the surface tessellations can adversely affect surface curvature and must therefore be carefully employed. Two case studies of sized based meshing applied to initial structured meshes elaborate on the above concepts. The mesher developed for this effort is discussed later but its results are used here for illustration purposes. For the first study, a quadrant of a truncated cone was generated with a structured mesh made of up quad facets that was then split into triangles. The resulting ?forced? unstructured mesh is shown in Figure-5.1.1. The initial sizing of the facets was set up so that the base of the cone had facets of appropriate edge-sizes. Note that this initial surface has all the inherent advantages and disadvantages of a mapped structured mesh. There are no surface perturbations and the feature lines are clean. Also note that the U-V mapping forces poor quality facets near the tip of the truncation edge of the cone and that the edge- sizes of the facets in this region are much smaller than the requested edge-size. Consider now the re-tessellated surface generated by the size-based mesher developed for this effort (Figure-5.1.1). The sizing distributions are much more uniform across the surface and especially near the tip of the truncation edge. However, note that 56 the feature lines on this surface are distorted and as a result there are surface perturbations (this can be evaluated in Figure-5.1.1 by the lighting distribution on the two surfaces from a single point source away from the surface). A comparison of the average edge-sizes and the final facet-count is shown in Table 5.1.1. The reduction in facet count is more than 50% while the improvement in edge-sizes created by the mesher is within 5% of user requested values compared with 40% on the original surface. Figure 5.1.1: Size based meshing case study for a truncated cone surface. Initial surface mesh (left) and re-tessellated mesh (right) Table 5.1.1: Size based meshing results for the truncated cone surface Initial surface Re-tessellated surface Number of facets 352 160 Average edge-size (expected = 0.2m) 0.28 0.21 57 A second case study was made on much more realistic application geometry. The Rutan Varieze civilian aircraft was run through the mesher and the results for the aircraft are shown in Figures-5.1.2(a) and Figure-5.1.2(b). In Figure-5.1.2(a), the original surface prior to re-tessellation is shown. This mesh was generated from a forced conversion of the underlying structured mesh which resulting in the starting surface of 108,040 facets. The geometry once again had clean feature lines and little to no surface perturbations. The large number of facets for this surface resulted from the very large needle-shaped facets distributed over the entire surface to allow effective pressure-gradient capture from a potential-flow pressure-solver. This surface was then run through the mesher with an aim to convert the surface into one more appropriate for a vorticity-solver. This meant that the sizing gradients could be improved in favor of reducing computation time and feature curves could be dispensed with where necessary. The result of this was process is shown in Figure- 5.1.2(b). As can be observed from this figure, the overall ?essence? of the geometry has been preserved by the mesher but at the cost of much more pronounced surface perturbations. Such a surface as shown in Figure-5.1.2(b) is not considered solver-ready for any pressure-field solver or any CFD solver but is usable by the current approach without loss of fidelity. The reduction in facet count is an indicator of reduced computation time outlined via Equation-3.3. The new facet-count is 10,036 facets (a reduction of over 90 %) compared with a structured mesh. The savings in computation time is over two orders of magnitude. 58 Figure 5.1.2(a): Initial pressure-solver mesh for the Varieze (108,040 facets) Figure 5.1.2(b): Re-tessellated vorticity-Solver mesh for the Varieze (10,036 facets) 59 A further illustration of the sizing distribution and gradients for this case study can be made via Figure-5.1.3. Here, the contour distribution of the mesh facet sizes is shown for the initial and final meshes on the Rutan Varieze. The comparison is qualitative but apparent. Figure 5.1.3: Mesh sizing distribution and gradients 5.2 MESH QUALITY DEFINITION The current effort quantifies mesh quality at the facet level and evaluates them at an area-averaged level for the overall surface mesh. At the facet level, the quality criterion (developed by the author) is defined as: ??????? = ????? ?? ?? ????? ??? ????? ?? ??? ???? ?? ??? (5.2.1) 60 The Circumcircle and Incircle for a triangular facet is readily defined in two dimensional Euclidean space based on the facet edge lengths and is thus applied to three- dimensional facets similarly. As a result of this, Equation-5.2.1 can be modified thus: ??????? = * ? ? ? 4? ( ? )( ? )( ? ) + ?( ? )( ? )( ? ) (5.2.2) Where, = (?12 + ?23 + ?31)2 (5.2.3) The limiting value for the quality ratio described in Equation-5.2.1 is 2.0 for an equilateral triangle. Note that this value represents the confluence of solver vortex ring requirements in that an equilateral triangle has the most appropriate placements of the centroid-based control points and therefore the resultant positive effect on solver convergence and reduced computation time. This lower reference value therefore serves as the metric for ?best-possible-quality? for a given facet on a surface. As a result, the mesher looks to this metric to determine which facets must be culled from a quality standpoint during the re-tessellation phase. That is, the mesher determines a range of quality above this reference value which it considers as recoverable and amenable for improvements. All facets having quality values beyond this threshold are culled from the mesh via edge collapsing (to retain conformity of neighbors). 61 Since the quality of a facet is geometrically defined, visual evaluation of overall mesh quality serves as a useful step for the current solver process. An example is shown in Figure-5.2.1 for the Lockheed-Martin F-35 geometry. Run through the mesher prior to re-tessellation, the ?quality-distribution? of the mesh can be visualized. In this way, mesh regions needing facet culling or re-tessellation can be identified and either manually repaired or run through the mesher. Figure 5.2.1: Unstructured mesh facet quality distributions 5.3 MESHING AUTO-REFINEMENT The automated unstructured mesher developed for this effort uses several invasive surface improvement techniques to render a size- and quality-controlled solver-ready mesh. The individual techniques applied at the facet or patch level are: a) Edge-collapsing b) Edge-Swapping c) Facet-splitting 62 Edge-collapsing is applied on those facet edges that are deemed too small to allow recovery of facet sizing and/or quality. An example is shown in Figure-5.5. The smallest edge in the facet is identified and the corresponding vertex connectivity map is determined from the global topology. The two vertex pairs are then collapsed and the topology adjusted accordingly as shown in Figure-5.3.1. Edge-collapsing serves to eliminate very poorly sized facets or ones with very poor quality (above the mesher recoverable threshold) while also helping to increase the overall size of the surrounding facets. Of the three surface modification techniques used in the mesher, it is the only one that increases the sizes of facets and reduces overall mesh size. Checks are also made to ensure that edges marked as feature edges based on their neighborhood connectivity are not subject to edge collapsing. This is done to ensure that no overall loss of geometrical ?essence? takes place. Figure 5.3.1: Edge collapsing on below-threshold facets. Before (left), after (right) Edge-swapping is a purely quality control technique and is applied iteratively in conjunction with the Edge-collapsing and Facet-splitting. The general principle is to 63 identify pairs of neighboring facets from the surface topology and identify their shared edge. If this edge can be flipped with resulting increase in quality of the two facets, the mesher allows this. Checks are made to evaluate the surface-normal angle across the shared edge and determine if this value is greater than the user defined threshold for active feature edge modifications. If so, the edge is not swapped and the feature edges are preserved. An example is shown in Figure-5.3.2 for two facet pairs. The resulting improvement in quality is visible. Note that the swapping subroutine actively looks to ensure that the final summed area of the two facets does not increase in order to prevent flipping edges that would lead to topologically inconsistent or self-intersecting facets in the post swapping phase. Note that edge-swapping is the only active quality-control technique in the mesher?s tool set. All other techniques may improve quality based on secondary criteria, but otherwise enforce no checks for quality. Figure 5.3.2: Edge swapping on facet pairs, before (left) and after (right) Facet splitting takes place in one of two ways: at the facet level or at the patch level. At the facet level, the mesher identifies facets that are above the allowable sizing 64 threshold and marks them for splitting. Topological sharing information for the three edges of the facet are recovered from the global map and the neighboring facets are also marked for splitting along the shared edge to enforce mesh conformity. Note that the splitting, whether at the facet level or the patch level, is designed to always enforce conformity on the surface unless none exists initially. Once this information is extracted, the mesher splits the edges in digital steps, i.e. the splits are dual along each edge and are not matched for exact sizing. This is enforced in order to prevent the formation of very poor quality, needle-shaped facets that can prove irrecoverable by the edge-swapping or edge-collapsing routines. Figure-5.3.3 shows an example of splitting applied by the mesher at the facet level for a simple cube designed with just 12 facets initially. Note the substantially higher facets in the resulting mesh in Figure-5.3.3. This is an inherent feature of facet splitting: there is always an increase in mesh size. As a result, only patches where the solver requires greater refinement are subjected to these routines. Note also that facet splitting does not respect feature edges in its effort to maintain mesh conformity. Figure 5.3.3: Individual facet-splitting, before (left) and after (right) 65 Facet splitting at the patch level is more advanced and sizing-driven. In this routine, the mesher actively looks for topological information and uses this to identify patches of the initial tessellation that can be suitably split along the feature lines. Once identified, it splits the feature lines in exact sizing steps (unlike the digital splits used in the facet-level splitting) so that feature-aligned surface meshes can be extracted. This is illustrated for an application made to the F-35 upper-fuselage surface from Figure-5.2.1 and is shown in Figure-5.3.4. This is the only technique in the mesher?s tool set to reduce overall feature-curve distortion in the final mesh compared to the initial surfaces. Figure 5.3.4: Patch-splitting, before (left) and after (right) Pierced facets are corrected by the application of splitting at facet level for the two intersecting pairs. The splitting point inside the facet is determined by the intersection of the edges of the opposite facet with the plane of the other facet. Once this point is evaluated, it is added to the global vertex map and the facet is split into three new facets about this point. This process retains the overall mesh conformity and the global edge-connectivity map is added upon but otherwise not changed. This process is applied 66 on both facets of the pair and the process is repeated until no further piercings exist on the pair. An example of this is shown in Figure-5.3.5. Here, a flat finely-tessellated surface is shown cutting a coarsely tessellated cube. This results in numerous piercings as visible from Figure-5.3.5. The post piercing-correction results are also shown and the increased facet count can be seen as well as the cutting edges that are not visible. The new facets are mainly of poor quality and further iterations of the solver are suitable candidates for edge-collapsing and edge-swapping. Figure 5.3.5: Piercing correction, before and after (cutting surface hidden) All of these techniques are applied in an iterative loop that spans every facet on the original mesh and changes are recorded. The mesher continues to operate on the mesh until it records no further changes applicable to the geometry at which points it terminates operations and moves the mesh to the solver for pre-simulation diagnostics. The liberal application of this mesher on the entire geometry can result in significant mesh-size increases, distortion of feature-edges and de-featuring of smaller components (that fall 67 entirely below the user specified sizing threshold). This can, of course, be used to advantage so that non solver-ready meshes with high geometrical fidelity can be reduced to a solver-ready state around the essence of the overall geometry. A case study of the auto-mesher is provided here for the F-16A geometry. Figure- 5.3.6 shows the initial geometry as run through the solver diagnostics. Facets marked in red are pierced and can cause instant solver divergence. Blue edges are feature edges evaluated by the mesher through topological probing. Note that feature edges have been detected on this geometry even over relatively smooth surfaces because of the underlying non-manifold edges that are not visible in Figure-5.3.6. Edges that are non-manifold are instantly marked as feature edges by the mesher and as such either need to be removed manually or the auto-mesher can leave them in place. Green edges are free and represent the perimeter of free-edges and non-conformal interfaces. Facets marked in orange show proximal facets that will also cause instant solver divergence. It is easily visible from Figure-5.3.6 that the initial mesh is highly articulated and contains a lot of electronics antennae and other empennages that are below the sizing threshold of the rest of the geometry and can be dispensed with (each of these empennages can cause solver divergence since they are not topologically consistent with the rest of the geometry; i.e. they pierce other facets or are in numerically bad proximity of each other). This repair can be accomplished via laborious manual means or through the use of the automated mesher developed here. The results are shown in Figure-5.3.7 and the quantitative mesh comparison is presented in Table 5.3.1. 68 Figure 5.3.6: F-16A initial geometry pre-solver topological probing results Figure 5.3.7: F-16A case study results for automated mesh refinements Table 5.3.1: F-16A automated mesh refinement results Initial Mesh Re-tessellated Mesh Number of facets 4504 13964 Average facet quality 8.3491 2.2253 Worst facet quality 4579.36 106.52 69 5.4 PATCH INTERFACES AND CONFORMITY For this effort, patch interface conformity is defined as the connectivity of facets at the feature edges of a surface topology. These interfaces are created during the CAD phase of the geometry design process but can also result from user defined feature curves as well as brought about within an MDO environment by the intersection of topologically independent bodies such as wing-fuselage intersections, wing-pylon intersections etc. Applied to potential-flow solvers at the surface level, they are the equivalent of CFD solver volume mesh partition interfaces and have philosophically similar effects on the solver convergence. These interfaces are handled in one of two ways: a) Mesh side handling b) Solver side handling In the mesh side handling, the mesher (volumetric for CFD and surface for potential-solvers) is responsible for generating conformal meshes along the interfaces. The solver simply inherits the connectivity information from the mesher during the solution phase. Maintaining conformity on the mesher side enforces mesh quality penalties that can adversely affect the solution. In the solver side handling, the interface conformity is ignored by the mesher and therefore meshes the topological patches in a standalone manner. The solver is then left to determine the cross-connected facets and neighborhood connectivity. This can result in increased computation time and increased solver oscillations depending on the quality of 70 the interfaces. The choice of either of these two approaches to interfaces affects the solver significantly. The vorticity-solvers and the Method of Integrated Circulation are largely insensitive to the surface presence or absence of mesh conformity at the patch interfaces with regard to the final aerodynamic loads. However, the solver convergence (and thus computation time) is still dependent on mesh conformity. The solver developed for this effort handles the patch interface conformity on the solver side. As such, the mesher is allowed to identify and mesh patches in standalone format. This reduces the complexity of the application of this approach to an MDO environment whilst improving flexibility in the design space since each patch can be meshed individually. Given the first-order approach of the vortex-ring solver applied to the facets, it was determined that the solver is generally stable around non-conformal patches under controlled sizing gradients (applied by the surface mesher described in previous Sections). A case study for the stability of the solver under conformal and non-conformal patches is presented here to illustrate the effects. The Soviet Mig-21 interceptor geometry was chosen for this study as shown in Figure-5.4.1. The wing surface was made conformal for one case and non-conformal for the second case (Figure-5.4.1) and both meshes were run through the current solver under identical conditions. The solution for vorticity and the relaxed wake strands were seen to converge in both cases and were found to be very close to each other (see Figure-5.4.1). The final solution therefore had 71 no apparent effect by the conformity of the wing-fuselage intersection. The integrated aerodynamic loads evaluated using the Method of Integrated Circulation seems to confirm this assumption in that far field effects are indeed diminutive on the final results. Figure 5.4.1: Mig-21 mesh conformity case study results However, such visual inspection of the vorticity distribution and the qualitative comparison of the integrated loads can be deceptive to the overall effect on the solver. To illustrate this, the convergence history for the above case study is provided in Figure 5.4.2. Plotted as a function of iterations, the dramatic effect of mesh conformity is easily visible on both the vorticity convergence as well as the wake-strand stabilization. From the figure it can be seen that the conformal design undergoes substantially more damped oscillations compared to the non-conformal mesh and the convergence takes place at a more rapid rate. 72 Figure 5.4.2: Mig-21 case study convergence history results Similar results have been obtained from other geometries as well that suggest that solver stability is adversely affected by the non-conformity of patch interfaces but otherwise have little effect on the overall fidelity of the integrated aerodynamic loads. 73 6. SOLVER EFFICIENCY ENHANCEMENT The basic vorticity-solver methodology utilizing structured meshes does not offer significant flexibility to improve its basic performance. However, with the current approach making it possible to apply unstructured meshes, several important techniques have been developed to take advantage of the increased solver and mesh flexibility to improve solver efficiency, stability and to reduce computation times. The novel techniques developed for this effort are: a) Fluid-Dynamic Mesh Repartitioning b) Persistent Geometrical Induction Matrix 6.1 FLUID-DYNAMIC MESH REPARTITIONING The unstructured surface meshes post CAD-driven tessellation and post auto- meshing (as described in Section-5 of this document) can render a diverse spread on the surface facet distributions with regard to numerical space. In other words, the unstructured nature of the mesh operations, whether part of CAD design, manual surface repair or the auto-mesher, leads to facets that are numerical neighbors being widely separated in physical space. This is illustrated in Figure-6.1.1 for clarity. The F-16A 74 geometry used here is illuminated for numerical neighbors using the contour coloring and shows the mesh post-CAD and post-auto-mesher. Figure 6.1.1: Rationale behind FDMR The widely spread nature of the numerical neighbors are readily appreciated from Figure-6.1.1. When passed to the solver in this condition, the iterative patches are scattered and have physical neighbors that are not part of the current iteration?s solution (and hence have vorticity values frozen from the previous iteration). As a result, the solver undergoes substantial stress. Once the solver stabilizes and the oscillations dampen out, the vorticity and wake solutions converge quickly to their final form. It has been observed from several geometries that the oscillation phase of the solver for a typical mesh as shown in Figure-6.1.1 constitutes about 40-50% of the total solver iteration time. If these oscillations could be damped out prior to start of solution, the resulting convergence would be substantially faster and more stable. Further, if the solver partitions could be organized so that partitions containing the leading edges and 75 sharp nose sections were evaluated first, the ?pass-on? effect of the solver during each iterative loop could lead to quicker convergence as well. To resolve this issue, the Fluid-Dynamic Mesh Repartitioning (FDMR) concept was developed and applied to the solver. This process repartitions the solver-ready mesh developed by the auto-mesher or provided by the user and reorganizes the numerical distribution of the mesh to match the free-stream direction. The sorting algorithm used here is a standard bubble-sort routine and the direction of the bubble-flow is in the free- stream direction. In some cases it is beneficial to sort the partitions in the direction normal to the flow and this has also been developed and added to the solver as a user defined option called Contra-FDMR. Figure 6.1.2 shows a case study for the application of the FDMR concept on the Lockheed-Martin F-35 geometry. This geometry was passed directly from the CAD to the solver post-FDMR (i.e. no auto-meshing refinement was used). As a result, it retains the original CAD tessellations during the solver run. This was done to ensure that the only change between various runs was the numerical re-ordering of the solver partitions. An example of FDMR application to the geometry is given in Figure-6.1.2 for flow in the vehicle X direction (long axis). The comparison of the solver partitions is clear and it is seen that the FDMR partitions provide the required ?pass-on? numerical behavior that is anticipated to give desired reductions in solver computation times. 76 Figure 6.1.2: The F-35 geometry before and after FDMR The geometry was run through the solver without FDMR, with FDMR and with Contra-FDMR for otherwise identical conditions. The vorticity solution for the geometry with FDMR is shown in Figure-6.1.3 and shows a clean convergence. Figure-6.1.3 shows that the FDMR solver convergence is an order of magnitude lower than the standard mesh and the Contra-FDMR results. In terms of computation time, the results are even more significant. The FDMR solution reaches the same convergence level of the non- FDMR solution (reached in 25 iterations) within 4-5 iterations. That is a reduction of ~80% in overall computation time during vorticity convergence phase of the solution. It should be noted that the loads evaluation phase for all three cases would remain unaffected by the application of FDMR. Also, FDMR imposes its own time penalty based on the mesh size. It is usually between 1-5% of the overall solver run time. Combining all of these factors, the effective reduction in solver time for the FDMR case is between 60- 70%. This value is still very significant and bears importance during MDO runs. 77 Figure 6.1.3: Convergence history for the F-35 solution (AOA = 5?, 1 m/sec flow) 6.2 PERSISTENT GLOBAL INDUCTION MATRIX As outlined in Section-4, the velocity induction coefficient matrix is identical during all iterations in the absence of geometrical warping or morphing. Under these conditions, it is possible to reduce the evaluation time for this matrix by storing it in the solver memory allocations after the first complete iteration. This approach reduces the computational expense for the solver and it merely has to access the predetermined memory locations to evaluate the geometrical induction coefficient of a given facet- 78 bound vortex ring on any other facet. This coefficient is then multiplied by its vorticity value from the latest iteration. Note that this does not apply to the relaxed wake strands since their spatial relaxation does not allow persistent storage of their induction effects of the geometry control points. The wake coefficients are evaluated incrementally despite the associated computational expense. The wake evaluations therefore become the limiting condition for the scalability of this feature for large geometries. 79 7. BOUNDARY PHYSICS The unstructured nature of the surface mesh used in the current effort makes the application of external boundary physics on local regions a simple process. Single facets can be marked individually or as part of topological patches. For the solver developed as part of the current effort, there are four main boundary conditions and associated external physical conditions that are applicable on the surface facets: a) Slip walls b) Velocity Inlets c) Trailing edges d) Symmetry 7.1 SLIP WALLS Slip walls are the default boundary conditions for all facets in the solver for the current effort as a carryover from all legacy potential-flow solvers. The current application takes the form of the Neumann condition 2, 14 on the control point of each facet (placed at the centroid of the facet). This condition states that the fluid flow normal to the direction of the facet at the control point must be zero. Unlike the other sister 80 condition for slip walls based on the Dirichlet condition 14, the Neumann condition remains applicable for both thick and thin surfaces and thereby helps maintain solver robustness 14. Legacy pressure-solvers were usually provided with the Dirichlet condition as the default slip-wall condition given the requirement for thick closed surfaces for pressure evaluations 8-11. Mixed-condition solvers were also used. However, vorticity solvers applying combination of thick and thin surfaces have always necessitated a Neumann condition approach 2. It also has a simpler implementation. The main disadvantage of the Neumann condition is the inability to carry or enforce any boundary layer information and thereby serves to decouple the viscous and inviscid effects, necessitating the requirement for a decoupled skin-friction load evaluation scheme, as outlined in Section-4. 7.2 VELOCITY INLETS Velocity inlets are simulated on the user-marked surface facets by relaxing the Neumann condition 14 on the control points. The relaxation of the Neumann condition allows the modeling of various types of inlet flows such as compressor fan faces of podded engines, air-intakes on the fuselage and cooling-intakes at various points on the geometry 14. The Neumann condition is provided the required velocity increment magnitude and direction vector components during its application in the solver. As a result, the solver ensures that no surface normal component of velocity other than the predetermined 81 external value can exist on the marked facets. It achieves this by modifying the local vorticity. Note that the Neumann condition only allows the modification and control of the surface normal velocity on the facet. It does not enforce any control on the local tangential components of the velocity (borrowing from its slip-walls lineage). As a result, relaxation of the Neumann condition to simulate velocity flow can only accurately portray the normal velocity components whilst ignoring the tangential components. An example of the application of velocity inlets is provided for the case of the Boeing F-18A Hornet (Figure-7.2.1). The geometry is modeled about the symmetry plane and has facets marked across the engine inlet next to the fuselage as shown in the figure. Two cases were run with the solver for this geometry for an AOA of 0?, symmetry plane applied and 25 iterations limit. The first case did not simulate the velocity inlet and treated the facets on the engine inlet as standard slip-walls (Figure-7.2.2). The second case relaxed the Neumann condition on the marked facets and simulated flow moving through the intakes (Figure-7.2.3). Notice the difference in the local vorticity distribution for the two cases. The under-wing vorticity is significantly affected by the presence or absence of the velocity-inlet and this translates to small changes in the integrated aerodynamic loads. The overall effect varies from geometry to geometry and must be exercised accordingly. 82 Figure 7.2.1: The F-18 geometry with velocity inlets and symmetry plane Figure 7.2.2: The F-18 geometry without Neumann relaxation Figure 7.2.3: The F-18 geometry with Neumann relaxation 83 7.3 TRAILING EDGES Facets at the trailing edges of wings, fins and stabilizers are marked accordingly so that the solver can determine the wake nodes on the mesh (as discussed in Section-4). As such, they remain the most important aspect of all potential-flow solvers, regardless of the solver methodology. The solver is set up to check marked facets for trailing edges as well as doing a global search on the surface mesh to determine trailing edge facets. In the former approach, the user marks the facets on the mesh manually to treat them as a trailing edge. This has the advantage of giving extensively customized solutions. However, for an MDO environment, this is clearly unsuited. Further, lacking the U-V mapping of structured meshes, the current solver cannot rely on pre-marked information for determining trailing edges on the mesh. As mentioned in Section-4, the primary advantage available from the relaxed strand model is the ability to evaluate wakes at the level of strands. As a result, the determining characteristics of the trailing edge facets can also be evaluated at the facet level in the mesh and then pooled together to determine the wake nodes. The current approach is based on the formulation of trailing edge metrics. The solver checks the facets in the entire mesh to determine which of the facets pass the testing criteria for trailing edges and automatically marks them as trailing facets. There are two main criteria for a facet to be a trailing edge: a) The facet must have at least one free edge or have a neighbor where the angle across the shared edge is greater than the wake definition threshold angle set 84 by the user (the stable value for this angle is 90? based on experimentations; Figure-7.3.1). b) The ?centrifugal-normal? for the edges that pass selection criteria (a) must make an angle with the free-stream velocity vector greater than the trailing- edge marker angle (the stable value for this angle is 75? based on experimentations). The centrifugal normal vector for an edge is the vector normal to the edge direction (along the winding of the facet) in the plane of the facet with a direction away from the facet centroid (Figure-7.3.2). Figure 7.3.1: Condition-1 for auto-detecting ?thick? trailing edges marked by the solver using topological feature angles 85 Figure 7.3.2: Centrifugal normal vectors for a general facet Facets that meet the above two conditions are processed for the wake generation criteria including the detection of possible wake nodes along the exposed trailing edge vertices (Figure-7.3.3). Figure 7.3.3: Typical results of trailing edge auto-detection. Thick wing example (left) and thin wing example (right) It should be noted that the solver is not restricted from marking more than one edge as a trailing edge. It is conceivable that up to two edges of a facet could simultaneously be trailing for a certain orientation and hence will both pass through the 86 auto-detection routines. The solver takes the existence of this possibility into account and lets the wake node detection routines determine whether a facet can have more than two wake nodes. In this aspect, the current approach attains a significant advantage compared to the pressure-solvers. The ability to detect thick and thin edges using the two selection criteria allows the solver to maintain robustness with both types of surfaces. This is yet another advantage of the current methodology. This allows the user to improve the efficiency of the solver by marking thick surfaces with very low thickness values as thin and eliminate and entire patch of facets from the solver without affecting the solver results. An example of such an approach is shown in Figure-7.3.4 for the Boeing F-18A Hornet. It has a geometry modeled from a combination of thick and thin surfaces whilst retaining overall geometrical (and hence the solution) fidelity. Figure 7.3.4: The F-18 geometry manifold and non-manifold surfaces 87 7.4 SYMMETRY PLANES For most standard applications, it is possible to identify a symmetry plane for the geometry with a normal vector perpendicular to the direction of the free-stream flow. This is especially true in those cases with non-morphing and non-rotating geometrical sections. Such cases can take advantage of the symmetry place boundary condition. The use of these planes significantly reduces the solver workload as well as memory requirements (in order to apply the persistent geometrical induction matrix feature discussed previously) by transferring over half of the mesh (and hence half of the wake) to the common induction matrix evaluated at the beginning of the solver iterations. When this condition is applied about the symmetry line, the facets retained are treated as though having a mirror image about the symmetry plane. Note that their winding direction does not change but the storage of vertices is modified accordingly. This allows the facet to induce its own effect spatially but also add to that the induction of its reflection on the other side of the symmetry plane. An example of the CPU savings from adopting this boundary condition is illustrated for the case of the DLR-F6 wing-fuselage geometry as shown in Figure-7.4.1. The original geometry has 11,038 facets. The effective facet count after symmetry plane has been applied is thus 5,519. Two case runs were done on the solver, with and without symmetry planes for 0? AOA and 1 m/sec free-stream flow. Both geometries were passed through the FDMR routines. The trailing-edges were detected automatically and the geometry was modeled entirely as thick. The resulting wake-field is shown in Figure- 88 7.4.2 and the results of the two runs are shown in Table 7.4.1 for the individual components that make up the overall solver CPU time. Figure 7.4.1: The DLR-F6 geometry with and without symmetry BC Figure 7.4.2: The DLR-F6 geometry with symmetry BC (AOA = 0?, 1 m/sec flow) It is easily established that the savings in CPU time for the overall run are substantial (64.46%) for the geometry run with the symmetry condition. Note also the reductions in FDMR time, global geometric induction matrix allocation and storage times 89 and the solver iterations time for the same number of iterations. There is reduction in all of these components, illustrating the advantage given by the symmetry plane condition. Table 7.4.1: DLR-F6 solver results comparison (AOA = 0?, 1 m/sec flow) No Symmetry BC Symmetry BC Number of facets 11,038 5,519 FDMR Wall Time (seconds) 2.58 1.02 Induction Matrix Evaluation Wall Time (seconds) 186.51 83.85 Solver iterations Wall Time (seconds) 454.14 143.70 Total Wall Time (seconds) 643.23 228.57 90 8. RESULTS This Section details the results of the solver runs made for the two main testing suites developed for this effort. The preliminary testing suite deals with simple geometries and shapes designed to test basic features of the solver. These include geometries such as spheres, flat wings etc. These geometries allow easy diagnosis of problem areas and help in identification of the limits of the current approach. They are also very useful to highlight the advances made by the current approach over legacy approaches and also to validate the new approach on well-defined theoretical cases. The advanced testing suite deals with real world applications and involves testing the solver under stressful mesh conditions (such as large mesh sizes, non-manifold surface intersections and complicated auto-mesher runs), with multi-physics external flows and testing limits of the physics of the potential flow philosophy. The geometries used here are real-world aerospace designs. Results have been collected for the aerodynamic load coefficients, solver run times, convergence histories and mesh statistics (wherever required). 91 8.1 PRELIMINARY TESTING SUITE The list of geometries tested in this group includes: a) Rectangle wings (various aspect ratios) b) Delta wings (various aspect ratios) c) 75? / 65? Double-delta wing d) 70? Arrow wing e) 70? Diamond wing f) 45? Sweep wing g) Unit sphere and other axisymmetric bodies 8.1.1 RECTANGLE WINGS Rectangle wings with aspect ratios of 3.0, 4.0 and 7.0 were designed and an unstructured mesh was generated for each of them. The geometry was designed as thin, free-surfaces. The results are shown in Figure-8.1.1 and Figure-8.1.2. The test results were compared with the original experimental results of Prandtl 23 from 1921 and are plotted alongside for comparison in Figure-8.1.1 and Figure-8.1.2. The results show generally good adherence to experiments with higher fidelity results for lower angles of attack. This is not unexpected since the new approach simplifies to the standard lifting line approach for rectangle wings in the extreme limiting conditions as disused in Section-4. At higher angles of attack, additional fluid physics such as flow separation begin to take place at higher angles and consequently affect the 92 lift and drag coefficients. The linear prediction curve then deviates from the experimental results as the experimental lift reduces while the experimental drag increases. Under these conditions the predictions overshoot the measured results. Figure 8.1.1: Lift coefficient versus angle of attack for rectangle wings 23 93 Figure 8.1.2: Lift versus Drag for rectangle wings 23 8.1.2 DELTA WINGS Unlike rectangle wings, delta shaped wings offer severe challenges for a vorticity- solver given the complicated physics even at modest angles of attack. Unless special care is taken to model the leading edge separation of vorticity as an additional relaxed wake proposed by Katz and Plotkin 14, or more advanced vortical sheet models are used such as 94 that done by Burkhalter and Purvis 24 for thin wing surfaces, the standard vorticity-solver is unable to capture the effects for delta wings except at low angles of attack 2. This approach enhances the application of such solvers to unstructured meshes and has devised load evaluation and wake generation techniques unique to such meshes, but otherwise remains true to the original philosophical concept. As a result, it is theoretically unexpected that the evaluation of loads over pure delta wings would be better. This is outlined in Figure-8.13 for delta-wings with various aspect ratios. The experimental data was evaluated for a Reynolds number of 7?105 by Schlichting and Truckenbrodt 25 in 1969. The experimental wings were modeled for thickness t=0.12c. It is easily seen that the experimental data is non-linear as a result of the physics involved. The current approach however can only match the linearized physics, and hence fits a linear lift-curve slope through the required data. The linearization deviates from experimental results at higher angles of attack and the error between prediction and experiments is significant within that region. 95 Figure 8.1.3: Lift coefficient versus angle of attack for delta wings 25 8.1.3 DOUBLE DELTA WING 75? / 65? In addition to the standard delta wings described in Section-8.1.2, several more wing shapes modified from the standard delta were also tested. One such model was the 75? / 65? double-delta wing. The plan-form of this geometry is shown in Figure-8.1.4 post auto-mesher. The nature of the unstructured mesh is also visible in this figure. 96 The results for this geometry were tested against the experimental results of Wentz and Kohlman 26, 29 from 1968. Additionally, the predicted data was also compared with the empirical equations based on the leading-edge suction analogy developed by Polhamus 27, 28 for high sweep angle wings. The comparisons are shown in Figure-8.1.5. It is seen that the behavior of the predicted data is similar to that of delta wings in that the current approach attempts to put a linearized prediction curve against non-linear data. The general fitting of the predictions is good for lower angles of attack similar to that for delta wings. Figure 8.1.4: Unstructured mesh on the 75?/65? double-delta 97 Figure 8.1.5: Lift coefficient versus angle-of-attack for the 75?/65? double-delta 26-28 8.1.4 ARROW WING 70? Another geometry modified from the original delta wings is the 70? delta wing as shown in Figure-8.1.6 post auto-mesher along with an overlay of its non-conformal unstructured mesh. The predictions for this wing are compared with the experimental results of Wentz and Kohlman 26 from 1968. Similar to Section-8.1.3, the Polhamus arrow empirical model 27-28 predictions have also been overlaid on the data for comparison. The arrow wing predictions are substantially poor compared with the double-delta and the delta wings even for the limited range of angles as seen in Figure- 8.1.6. The predictions are comparable to the experimental results for angles up to 4?. The substantial non-linearity of the arrow wing even at lower angles renders any such 98 potential flow approach ineffective. Burkhalter and Purvis have managed to get superior results compared to the current approach for such wing geometries using continually distributed vortical sheets. Figure 8.1.6: Unstructured mesh on the 70? arrow wing 99 Figure 8.1.7: Lift coefficient versus angle-of-attack for the 70? arrow 26-28 8.1.5 DIAMOND WING 70? The final modification to the delta wings was in the form of a 70? diamond wing as shown in Figure-8.1.8 post auto-mesher. The experimental data used for comparison for this geometry also follows from the work of Wentz and Kohlman 26 from 1968. The Polhamus equations 27-28 modified for a diamond wing are also used for comparison. The results are shown in Figure-8.1.9. The general comparison is good for the range of angles shown with the predictions underestimating the measured data. This is a trend that has now been established as common for all the delta-wing variations seen previously and in contrast to those for the rectangle wings, which always overestimated measured data at higher angles of attack. 100 Figure 8.1.8: Unstructured mesh on the 70? diamond wing Figure 8.1.9: Lift coefficient versus angle-of-attack for the 70? diamond wing 101 8.1.6 SWEEP WING 45? In contrast to high taper-ratio wings such as the delta, diamond and arrow wings, the current approach predicts very good results for low taper-ratio geometries and does so for other geometries that have physics that are captured by vorticity-solvers. This was seen in the limiting case for a rectangle wing but is not limited to such basic geometries. The case presented here for a 45? sweep wing is another example where the current approach predicts the aerodynamic performance in close accordance with measured data. The geometry used here is a standard zero taper-ratio rectangle wing with an aspect ratio of 5.0 and 45? sweep on the leading and trailing edges. The experimental data used for comparison here is obtained from the work of Weber and Bremner 30 in 1958 and is shown in Figure-8.1.10 and Figure-8.1.11. Figure-8.1.10 shows the very close approximation of the predictions alongside the measured lift coefficient for the wing for the entire range of measured angles of attack. Similarly, the sectional lift coefficient plotted as a function of span-wise location also shows remarkable estimations alongside the measured data. Note that the sectional lift distribution is nonlinear for such wings 30, 31 but the current approach is able to match that nonlinearity accurately. Some noise is seen inside the experimental data as well as the characteristic numerical noise in the sectional circulation data of the current approach (which manifests itself as the sectional lift coefficient using the Kutta-Joukowski formulation). 102 Figure 8.1.10: Lift coefficient versus angle-of-attack for the 45? sweep wing 103 Figure 8.1.11: Normalized sectional lift coefficient for the 45? sweep wing 8.1.7 UNIT SPHERE AND AXISYMMETRIC BODIES In addition to the standard wing geometries, axisymmetric bodies were also tested for the classical potential flow results using the current solver. Two of these bodies are presented here: a sphere and an ogive body. The sphere was designed with a unit radius and meshed unstructured as shown in Figure-8.1.12. The theoretical lift coefficient from such a body without wake description is zero. The geometry was tested from 0? to 90? angle of attack and the results are shown in Figure 8.1.13. The overall prediction is found to be well within acceptable limits around the theoretical prediction. There is some 104 numerical error obtained from both the solver as well as the facetization of the analytical sphere surface and its effects are random as seen in Figure-8.1.13. Figure 8.1.12: The unit sphere Figure 8.1.13: Lift coefficient versus angle of attack for the sphere 105 The ogive body used for the tests is shown in Figure-8.1.14. This body was generated with a center body radius of 0.5 meters and an overall length of 4.0 meters. The unstructured mesh generated on this body shows some of the limitations of the unstructured mesh discussed in previous Sections. The nose of the ogive is poorly discretized by the mesher but has very good quality facets. This is helpful to the solver but poor for overall mesh fidelity. Using an structured mesh for such a body would have increased the mesh fidelity by retaining the appropriate feature curves but would have resulted in solver oscillations and higher numerical errors the loads. As such, the predicted results for the lift coefficient are much better than the sphere as seen in Figure- 8.1.15. The numerical error is significantly lower and overall drives the point that the current approach requires only the ?essence? of the underlying geometry for it to provide accurate results. This will prove very useful in the following Sections that deal with the advanced geometries. Figure 8.1.14: The ogive axisymmetric body 106 Figure 8.1.15: Lift coefficient versus angle of attack for the ogive body 8.2 ADVANCED TESTING SUITE The list of geometries tested in this group includes: a) General-Dynamics F-16A b) McDonnell Douglas F-18A c) DLR-F6 commercial transport aircraft d) NASA Canard Fighter e) Rutan Varieze 107 8.2.1 GENERAL-DYNAMICS F-16A The General-Dynamics F-16A Fighting Falcon was originally designed as a multipurpose strike aircraft in the lightweight fighter category in the 1970s 32. Its original design objective was superior subsonic-cruise lift-to-drag ratio to allow effective mission radii with a variety of air-to-ground and air-to-air weapons. The design (Figure-8.2.1) has leading edge extensions for the main wing to provide controlled vortex lift and increase lift at high angles of attack 32. The single vertical stabilizer was chosen over the two surface options unlike the F-15 to reduce buffeting from strake vortices at high AOA. The engine intake is place below the nose to avoid gun gas ingestion 32. The overall wing-body intersection is blended for superior aerodynamic performance and increased internal volume for fuel (and hence range). Figure 8.2.1: The Lockheed-Martin F-16 32 This design was tested using the current approach by the conversion of a flight- simulation geometry used by NASA into a solver ready mesh. The original geometry is 108 shown in Figure-8.2.2. This represented the starting point in the conversion process for the solver. Figure 8.2.2: Original geometry for the F-16A post CAD design This geometry was then cleaned up using manual repair tools to remove geometrical features that would not contribute to the solver fidelity but would otherwise increase the overall mesh size and induce solver instability and increase overall computation time. The features removed in the manual repair process included the weapons stores, the under-wing pylons, the gun ports and the electronics antennae. The engine intake and exhaust were closed off with additional facets and the internal geometry removed. Further, the thin nature of the horizontal stabilizer allowed for the collapsing of the upper and lower surfaces into a mean camber line thin surface. 109 Following the geometrical de-featuring, the model was passed through the auto- mesher to improve overall mesh quality and sizing gradients. The need from this can be seen from the auto-mesher results shown in Figure-8.2.3. Additional feature curves were marked manually to ensure the auto-mesher did not de-feature the surface curvature in the wing-body intersection regions as well as the cockpit-fuselage intersection regions (this allowed the preservation of the classic bubble-shaped canopy lines of the F-16). Note from Figure-8.2.3 that the auto-mesher induced surface perturbations at several points on the geometry. These were manually cleaned after the meshing process. The user of the auto-mesher increased facet count substantially as seen from the statistics in Table-8.2.1. However, the overall mesh quality was much improved. This was considered a good tradeoff between computation time and solution fidelity. The boundary physics for the geometry included slip walls at every point on the geometry except for the intake and exhaust facets that were marked as velocity inlets and outlets. The trailing edge detection was manual and the facets were marked prior to solver initialization. The geometry was evaluated using a symmetry plane. 110 Figure 8.2.3: Auto-mesher results on the F-16A Table 8.2.1: Auto-mesher results for the F-16A Initial Geometry Refined Mesh Number of facets 4504 13964 Average facet quality 8.3491 2.2253 Worst facet quality 4579.36 106.52 111 The solver was run for a 10-4 convergence setting and the Trefftz plane was marked 6 meters behind the horizontal stabilizer trailing edges. This was considered a safe point for the termination of the wake as a tradeoff between wake fidelity and computation time. The results of the vorticity-solver for the geometry are shown in Figure-8.2.4 for an AOA of 10?. The wrap-around of the main-wing wake strands near the fuselage is easily visible. Figure 8.2.4: Vorticity distribution on the F-16A The aerodynamic loads were evaluated using 32 ICL planes (using the average edge-size of the mesh) as shown in Figure-8.2.5. The sectional lift-coefficient as a function of span-wise location for an AOA of 10? is shown in Figure-8.2.6. The loss of lift caused by the presence of the fuselage is easily visible as well as the additive effects of the horizontal stabilizer. 112 Figure 8.2.5: ICL distribution about the symmetry plane on the F-16A Figure 8.2.6: Lift coefficient span-wise distribution about the symmetry plane The integrated lift data was compared with experimental results obtained by Nguyen et.al 32, in 1979 as part of a NASA effort. The flight conditions were highly 113 compressible at Mach 0.9 and the predicted results were corrected for compressibility using the correction terms described in Section-4. The final results are shown in Figure- 8.2.7. The overall comparison is found to be very favorable along the measured data. There is some deviation at the lower and upper spectrum of the data as a result of slightly lower slope of the predicted lift-curve. It is reasoned that this is the result of inaccuracies in the source geometry. Figure 8.2.7: Lift coefficient versus angle of attack for the F-16A 32 The integrated drag predictions were compared with measured data obtained from the work of Webb and Kent 33 in 1977 and were also collected for Mach 0.9 in flight. The predicted data fits well within the flight-test data but the induced drag is seen to be 114 underestimated at higher angles of attack. The effects of compressibility on the predicted data could also be a factor in this discrepancy. Figure 8.2.8: Lift versus drag for the F-16A 33 8.2.2 MCDONNELL-DOUGLAS F-18A The second advanced geometry used in this testing suite was the McDonnell- Douglas F-18A (Figure-8.2.9). The pre-solver processing applied on this geometry was very similar to that adopted for the F-16 including the manual de-featuring of extraneous external and internal geometry objects prior to auto-meshing. The final mesh is shown in Figure-8.2.10. The geometry was modeled with 7,857 facets about the symmetry plane. 115 Given the nature of the experimental results available for comparison, the engine intakes and exhausts were not modeled with a velocity inlet boundary. The horizontal and vertical stabilizers were converted to thin surfaces but the main wing was modeled as a thick body merged with the similarly thick LERXs (visible in Figure-8.2.10). Figure 8.2.9: The McDonnell Douglas F-18 Figure 8.2.10: Auto-mesher results for the F-18A The solver was run for a 10-4 convergence setting and the Trefftz plane was marked 5.0 meters behind the horizontal stabilizer trailing edges. This was considered a 116 safe point for the termination of the wake as a tradeoff between wake fidelity and computation time. The results of the vorticity-solver for the geometry are shown in Figure-8.2.11 for an AOA of 1?. Figure 8.2.11: Vorticity distribution on the F-18A (vorticity units = m2/sec; AOA = 1? and free-steam velocity = 1 m/sec) The aerodynamic loads were evaluated using 44 ICL planes (using the average edge-size of the mesh) as shown in Figure-8.2.12. The sectional lift-coefficient as a function of span-wise location for an AOA of 1? is shown in Figure-8.2.13. The figure illustrates the important additive and subtractive effects between the main wing, the 117 stabilizers and the fuselage. Moving from the wingtip inwards, the additive effect of the horizontal stabilizer is seen initially, spiking the overall sectional lift. This is then subtracted from by the tilted vertical stabilizer which caused the downward spike. The LERXs then add to the lift and the curve spikes up again and finally drop off above the fuselage, which does not produce any useful lift at such low angles. Figure 8.2.12: ICL distribution about the symmetry plane on the F-18A 118 Figure 8.2.13: Lift coefficient span-wise distribution about the symmetry plane The integrated lift data was compared with experimental results obtained by Banks 34 in 1988. The experimental conditions were slightly compressible at Mach 0.33 and a Reynolds number of 5.33?106. The predicted results were corrected for compressibility using the correction terms described in Section-4. At this point it is perhaps useful to compare the results of the experimental results of Banks with the vorticity solver used by Browne and Katz 35 in 1990 as well as the current solver developed for this effort. The geometry used by Browne and Katz is shown in Figure-8.2.14(a). The fidelity of the geometry is quite poor compared with geometry used for the current solver (Figure-8.2.14(b)). This is understandable given the computational hardware limitations in 1990. However, moving past the geometrical fidelity, several philosophical differences can be outlined between the past work and the current methodology. Primarily, the model developed by Browne and Katz 35 is a structured mesh with quad valency panels. This is in line with the approach taken by 119 legacy Vorticity solvers. For geometries as advanced as the F-18, this imposes a penalty with the solution fidelity compared with the unstructured mesh used for the current approach. Figure 8.2.14(a): Panel model of the McDonnell Douglas F-18A 35 with 1,336 panels (effective 2,372 facets for current solver) Figure 8.2.14(b): Unstructured model of the McDonnell Douglas F-18A with 7,853 facets used in the current solver This penalty between the two solver approaches is shown in the comparison of the results next to the NASA Ames wind-tunnel results by Browne and Katz 35 in 1990 (Figure-8.2.15 and Figure-8.2.16). As seen from the lift data, the comparison of both solvers is favorable at lower angles of attack (AOA < 12?). However, the vorticity-solver 120 used by Browne and Katz 35 under-predicts the induced effects while the current approach over-predicts the induced loads beyond AOA=15?. The drag results amplify this discrepancy substantially in favor of the current approach. It is seen from Figure- 8.2.16 that the current approach predicts the drag results much better than the older solvers at higher AOA. Even so, beyond AOA=20?, the current approach also begins to under predict the drag as nonlinear physics become dominant on the upper surfaces of the F-18 LERXs. Figure 8.2.15: Lift coefficient versus angle of attack for the F-18A 34-35 121 Figure 8.2.16: Lift versus drag for the F-18A 34-35 8.2.3 DLR-F6 COMMERICAL TRANSPORT AIRCRAFT The third geometry used in this testing suite was the AIAA DLR-F6 (Figure- 8.2.17) 36-38. The pre-solver processing applied on this geometry differed from the previous cases in that the original CAD tessellations did not have to be de-featured or any components removed other than the volumetric bounding box. This is because the DLR- F6 geometry is part of the AIAA drag prediction workshop 36-38 testing suite and as such is designed in CAD for use in CFD. Unlike the F-16 and F-18 geometries, which were designed as visualization meshes for simulators, this geometry is a converted CFD 36 model used in the current solver. 122 Figure 8.2.17: The DLR-F6 geometry on a mount in the DLR wind-tunnel 36 The initial CAD tessellations are shown in Figure-8.2.18. The geometry is already in an unstructured state and is used as such. The initial mesh was topologically consistent and the bounding box mesh about the symmetry plane was split accordingly 36 as seen in Figure-8.2.18. This was subsequently removed for testing in this approach since volume meshes were not needed. The presence of underlying CAD also allowed for a unique test of the CAD projection abilities of the auto-mesher which resulted in a good quality, solver-ready surface mesh as shown in Figure-8.2.19. This final mesh had 11,038 facets about the symmetry plane. On the DLR-F6, the podded engine is hollow inside and follows the model design for the wind tunnel 37-38 (Figure-8.2.17). This negated the need for velocity inlets to be applied on the engine intake and exhaust. Further, the geometry is simplified and there are no horizontal and vertical stabilizers on the geometry. No surfaces had to be collapsed to remove their thicknesses. True to the original CFD geometry, the model was tested with a symmetry plane boundary. 123 Figure 8.2.18: Initial AIAA CAD tessellations for the DLR-F6 36 Figure 8.2.19: Auto-mesher results for the DLR-F6 The solver was run for a 10-4 convergence setting and the Trefftz plane was marked 1.0 meters behind the horizontal stabilizer trailing edges (the model overall length is 1.188 meters using the scaling applied on the wind-tunnel geometry). As such the Trefftz location corresponded to roughly one body length past the last facet on the mesh and this was considered a safe point for the termination of the wake as a tradeoff 124 between wake fidelity and computation time. The results of the vorticity-solver for the geometry are shown in Figure-8.2.20 for an AOA of 1?. Figure 8.2.20: Vorticity distribution on the DLR-F6 (vorticity units = m2/sec; AOA = 1? and free-steam velocity = 1 m/sec) The aerodynamic loads were evaluated using 46 ICL planes (using the average edge-size of the mesh) as shown in Figure-8.2.21. The sectional lift-coefficient as a function of span-wise location for an AOA of 2? is shown in Figure-8.2.22. As for the case of the F-18, the figure illustrates the important additive and subtractive effects 125 between the main wing and the fuselage. However, unlike the F-18 geometry, the much larger numerical noise is seen in the ICL vorticity distributions as a result of the low- density mesh on the wing, which results in much larger surface perturbations. These can be improved with further local refinements of the mesh but as is outlined later, the effect on the results is expected to be minimal. Figure 8.2.21: ICL distribution about the symmetry plane on the DLR-F6 Figure 8.2.22: Lift coefficient span-wise distribution about the symmetry plane 126 The integrated lift and drag data was compared with experimental results obtained in the ONERA 36 tunnels (in 2001) as well as with two Navier-stokes solvers (CFL-3D and OVERFLOW 37-38) using two ?medium? meshes and Overset mesh solvers 36-38. The choice of using medium meshes as a comparison point with CFD is deliberate. It is well understood that the above solvers working with ?fine? meshes result in superior results than the relatively coarser ?medium? meshes. However, it was discovered that the tessellations used in the current approach are similar to the ?coarse? meshes used in CFD solvers, even though from potential-flow standpoint, they represent very refined meshes. DLR-F6 results with similar coarse CFD meshes could not be obtained 36. However, medium mesh results are available 36-38 and used here to illustrate the philosophical differences in requirement between the two approaches using the same underlying geometry. The experimental conditions as well as those used in the CFD solvers were highly compressible at Mach 0.75 37. In both these cases, the wing pressure distribution showed supercritical performance and normal shock spikes were measured on the surfaces during experiment 36. This remains out of bounds of the limit of the potential flow solvers. However, using the precedence from the F-16A tests described previously, the tests on this geometry were continued after including the corrections for compressibility as described in Section-4. The overall comparison is found to be very favorable alongside the measured data as seen in Figure-8.2.23 and Figure-8.2.24. The lift data shows the Navier-Stokes solvers with medium meshes overestimating the measured data while the current solver is able to 127 match closer. The perturbations in the evaluated data for the current approach are due to the perturbations noted in the ICL vorticity distributions in Figure-8.2.22 and are entirely numerical errors. The drag data in Figure-8.2.24 shows similar trends as the lift data for the two Navier-Stokes solvers 36-38. The hypothesis by the authors is that there is an over- prediction in the surface pressures as a result of the coarser mesh. Since the aerodynamic loads are evaluated by pressure-summations, both the lift and drag are seen to be over- predicted. The drag predictions by the current approach shows a little discrepancy once again brought out by the vorticity perturbations in the ICL (drag data between 0.025 and 0.03 is perturbed adversely). Given the outer bounds of the physics involved, the current approach manages to get relatively close matches with the measured data. Figure 8.2.23: Lift coefficient versus angle of attack for the DLR-F6 36-38 128 Figure 8.2.24: Lift versus drag for the DLR-F6 36-38 8.2.4 NASA CANARD FIGHTER The final two geometries in this testing suite deal with stress testing of the relaxed wake-strand model. Both geometries have canards and main wings which involve wake interference in the loads as well as complex wake relaxation. The similarity between the geometries ends there, however. The NASA Canard Fighter 39 was designed in 1983 purely as a testing model for improving the understanding of such canard-wing interactions. It has a configuration of a high subsonic combat aircraft similar to the Swedish Viggen and in recent years the French Rafale and the European Eurofighter designs 39. Unlike these real world geometries, however, the NASA design was a 129 simplified model with straight rectangular section and primary focus was given only to the canard-wing aerodynamics 39, 41-43. The model of the Canard-Fighter (shown in Figure-8.2.25 as mounted in the Ames wind tunnels) used for this effort was reconstructed from the detailed dimensions and scaling information provided in the work by Stoll and Koenig 39 at NASA Ames in 1983. This reconstructed geometry is shown in Figure-8.2.26 following auto-meshing cleanup. The mesh had 2,383 facets about the symmetry plane and was specifically modeled to include thin main wings and canard. This is in accordance with the very low thickness of the NACA 65A004 airfoil 40 used on the NASA model for both the surfaces. The savings in computation time as a result of this action is significant (a reduction of 23% was observed). And as outlined later, the effect on the final results is not noticeable. Also in accordance with the wind tunnel model 39, the engine intakes and the base of the model were modeled with facets but otherwise were not included with a velocity inlet boundary. The wing and canard mesh patch interface with the fuselage is non- conformal and was not fixed. The solver adapted to this condition without any oscillations given the general alignment of the facets in those regions. 130 Figure 8.2.25: The Canard Fighter model on a mount at NASA Ames tunnel 39 Figure 8.2.26: Auto-mesher results for the Canard Fighter post CAD design The solver was run for a 1e-4 convergence setting and the Trefftz plane was marked 5.0 meters behind the main wing trailing edges (the model overall length as used in the solver is 5.053 meters). As such the Trefftz location corresponded to roughly one body length past the last facet on the mesh and this was considered a safe point for the termination of the wake as a tradeoff between wake fidelity and computation time. The 131 results of the vorticity-solver for the geometry are shown in Figure-8.2.27 for an AOA of 10?. Figure-8.2.27 is also a good illustration of the stability of the wake strands. For this angle of attack, the wake from the canard moves above the main wing by a significant distance. This is partly due to the higher elevation of the canard relative to the wing and partly because of the fluid flow angles. The high incidence of the flow also creates powerful tip vortices for both the canard and the main wing and is easily seen in Figure-8.2.27. The vortices leaving the canard are under the influence of the main wing and the fuselage and as a result they do not spiral in a manner similar to the main wing strands. The strand emanating from the canards near the root are also pulled above and around the top surfaces of the fuselage as a result of the local vorticity gradients caused by the sharp angles on the geometry. The aerodynamic loads were evaluated using 18 ICL planes (given the relatively coarse mesh) as shown in Figure-8.2.28. The sectional lift-coefficient as a function of span-wise location for an AOA of 10? is shown in Figure-8.2.29. As for the case of the F- 18, the figure illustrates the important additive and subtractive effects between the main wing, canard and the fuselage. Given the simplified geometry involved, there is very little numerical noise in the ICL distribution. 132 Figure 8.2.27: Vorticity distribution on the Canard Fighter (vorticity units = m2/sec; AOA = 10? and free-steam velocity = 1 m/sec) 133 Figure 8.2.28: ICL distribution about the symmetry plane Figure 8.2.29: Lift coefficient span-wise distribution about the symmetry plane The integrated aerodynamic loads were compared with the experimental results obtained from the work of Stoll and Koenig 39, 41-43 at NASA Ames in 1983. Further comparison is made with the nonlinear Vortex Lattice Solver (VLS) used at Ames 44 by the authors, as well as with the Boeing PANAIR pressure-solver also used to compare 134 with the experimental results 45 (Figure-8.2.30). The VLS solution simplified the geometry by removing the fuselage and extending the wing and the canard to the centerline of the geometry 44. It also maintained the surfaces as thin. The PANAIR solution modeled the geometry entirely as a closed and manifold body for a pressure- integrated solution 45. As seen from Figure-8.2.30 for the lift evaluations, the current approach closely approximates the wind-tunnel results. The VLS predictions are also close to the measured values. However, the PANAIR results deviate quickly from the measured data with increasing angle of attack. The PANAIR deviation is expected 45 given the low facet count and poor geometric fidelity of the model: as outlined earlier in this document, pressure-solvers remain limited by the fidelity of the mesh used. However, the reason the non-linear VLS solution matches experimental data so accurate is because it models the rollup of free vortices similar to the current approach?s relaxed strand model 44. More importantly, the model has the ability to simulate vortex separation from leading, side and trailing edges using free-separation vortices composed of discrete vortex segments 44. The current approach does not model free-vortex separation from the leading edges and this remains one of the major future objectives. The accuracy of the Ames non-linear VLS solver is a classic example of enhancing the simplified VLS solvers for improved accuracy and physics 14, 44. Figure-8.2.31 illustrates the results for the model with and without the canards on the geometry. The comparisons here are restricted only to the wind-tunnel results given the absence of similar data from the NASA references cited earlier. The current approach 135 accurately captures the presence and absence of the canard from low AOA up to the point of vortex breakdown above the main wing surfaces. Figure 8.2.30: Lift coefficient versus angle of attack for the Canard Fighter 44-45 136 Figure 8.2.31: Lift versus angle of attack for the Canard Fighter 44-45 8.2.5 RUTAN VARIEZE The final geometry used in this testing suite increases the stress on the relaxed wake-strand model by combining a canard-wing with a realistic fuselage and more complex surface blending. The model chosen for this test is the Rutan Varieze 46 which is a small civilian pusher design as shown in Figure-8.2.32. The testing goals remain the same as those for the NASA Canard-Fighter 39. 137 Figure 8.2.32: The Rutan Varieze The geometry is designed using the detailed CAD information provided by the work of Long 47 in 1985 as part of a NASA research effort. The CAD was modified to remove the undercarriage fairings and booms as well as the propellers. The propeller nose cone was smoothened and blended into the rear fuselage. The vertical stabilizers at the wingtip were made into thin surfaces. Smaller fairings and antennae were removed as part of the standard auto-mesher de-featuring. The final resulting mesh is shown in Figure-8.2.33. The wing meshes had significant perturbations but were found to be irrelevant to the overall solver fidelity, in line with the robustness expectations following the DLR-F6 geometry testing done previously. Unlike the Canard-Fighter, the vertical stabilizers were modeled as lifting surfaces given their location on the main wing tips 46. The under-fuselage engine air intake was not modeled as a velocity inlet and was closed off as a solid surface. The geometry was modeled with a symmetry plane boundary. 138 Figure 8.2.33: Auto-mesher results for the Varieze post CAD design The solver was run for a 1e-4 convergence setting and the Trefftz plane was marked 5.0 meters behind the main wing-trailing edges. The results of the vorticity-solver for the geometry are shown in Figure-8.2.34 for an AOA of 10?. As discussed for the same geometry in Section-4, the wake strands are seen to show complex behavior over the wing surfaces and next to the fuselage and require more iterations (and hence more computation time) than for other relatively simpler cases such as the DLR-F6 geometry. 139 Figure 8.2.34: Vorticity distribution on the Varieze (vorticity units = m2/sec; AOA = 10? and free-steam velocity = 1 m/sec) The aerodynamic loads were evaluated using 48 ICL planes (using the average mesh edge-size as discretization length) as shown in Figure-8.2.35. The sectional lift- coefficient as a function of span-wise location for an AOA of 10? is shown in Figure- 140 8.2.36. As for the case of the Canard-Fighter, the figure illustrates the important additive and subtractive effects between the main wing, canard and the fuselage. Unlike the Canard-Fighter, however, the Varieze has substantially more numerical noise in the ICL distribution and has convergence history behavior similar to that of the DLR-F6, which also had substantial ICL distribution scatter due to mesh perturbations. Figure 8.2.35: ICL distribution about the symmetry plane 141 Figure 8.2.36: Lift coefficient span-wise distribution about the symmetry plane The integrated aerodynamic loads were compared with the experimental results obtained from the work of Long at NASA Langley 47 in 1985. The flow is incompressible with only Mach 0.1 flow. The tests on the Varieze shed important light on the coupling effects of the relaxed wake flow over lifting surfaces and their corresponding effect on the overall lift generation 46-47. Tests were done on the geometry with and without canards and with and without main-wing 47. The results of these tests for the lift generated are presented in Figure-8.2.37. From Figure-8.2.37, it is possible to see the coupling effects clearly. The tests done with individual components (canard or wing) are plotted and in both cases the current approach matches them accurately within the appropriate AOA range prior to flow separation. However, a plot has been made to show the decoupled sum of the wing and canard lift coefficients. This curve is parallel to the experimental data for the coupled sum of the canard and the wing. This is the coupling effect of the canard wake moving 142 over the wing surfaces that causes the overall lift coefficient to reduce. When run with the complete geometry setup, the current approach captures the effect clearly and the matches the measured data closely. This test underlines the importance of (and need for) relaxed wakes for complex aerodynamically-coupled geometries such as the Varieze. Figure 8.2.37: Lift coefficient versus angle of attack for the Varieze 47 9. CONCLUSIONS A method for the application of vorticity based potential-flow solvers to unstructured surface meshes has been created through the use of a novel method for aerodynamic loads. This method is designed to maintain the advantages of the vorticity based solvers while attempting to remove the limitations of geometrical inputs associated with the philosophical approach. A discussion on the necessity and advantages of this approach has been presented. It was shown for the case of simple geometric objects such as the sphere that 65% reductions in solver iterations could result from the use of an unstructured mesh (with inherent higher mesh quality) versus a converted structured mesh. In addition, a reduction of ~54% in mesh size was obtained as a result of higher mesh quality. Overall, the combination of reduced mesh size and higher quality led to a reduction of about ~88% in overall solver CPU time. The theoretical underpinnings of the current methodology have been elaborated on in this document. A classical test case is presented in the thin elliptical wing for which an analytical solution is known from lifting line theory. Comparisons of the current numerical approach with this analytical solution show close approximation (error variation within 2.5-3%). The linearized nature of the underlying theoretical foundation for the current theory was shown using the F-18 geometry. In accordance with the stated 144 objectives of this effort, stability metrics were developed for the surface vorticity and the relaxed wake strands. These parameters have since been used to demonstrate the stability of the solver under stressful geometry conditions. To further capitalize on the improved possibilities offered by the use of unstructured surface meshes, a modification to the evaluation of skin-friction coefficients using surface vorticity was developed to bypass the limitation of existing approaches. The heuristic arguments in support of the current formulations have been presented and fit well within experimentally observed behavior. A novel unstructured wake-strand model has been developed to allow handling of wakes emanating from unstructured surface meshes and test cases show the stability of the vorticity shedding as well as the applicability of its novel features such as adaptive time stepping and vorticity decay models used to eliminate mathematical singularities. An attempt has been made to extend Potential-Flow solvers to the current industry surface meshing and numerical solver standards. This automated design pipeline is considered robust enough to allow integration with optimizers within an MDO environment no longer restricted in their geometrical design spaces. In accordance with the objective of this effort, metrics were developed for surface mesh quality and automated meshing tools were developed to convert input geometries into solver ready meshes. The effect of size based meshing and mesh quality on solution fidelity was tested. It was determined that surface meshes with uniform sizes and quality offered superior results and quicker convergence. RMS error reduction between 10-20% between higher and lower quality meshes. 145 Several improvements in solver methodology applicable in an unstructured mesh solver have also been developed. These included the FDMR and persistent induction matrix concepts. The use of the FDMR on the F-35 fighter geometry resulted in 80% reduction in solver iterations for the same convergence limits and an order of magnitude improvement in solver convergence (10-3 versus 10-2 for the non FDMR case) for the same number of iterations. Testing using the persistent induction matrix concept allowed reductions between 65-85% in CPU time for the same number of iterations for large surface meshes (>10,000 facets). Improvements to the boundary physics models have been made to take advantage of the unstructured meshes. Relaxation of the Neumann condition on user-specified facets or topological patches has made it possible to create velocity inlets. Additionally, the detection of the trailing edges has been automated, thereby removing a possible limitation compared to structured mesh solvers. Use of the Mirror-plane boundary condition has allowed substantial reductions in solver time for steady-state, symmetric flow geometries. Tests on the DLR-F6 geometry with and without the mirror plane boundary showed that reductions of up to ~65% could be obtained through the use of the mirror plane. Other geometries allow for varying savings in CPU time with the general trend towards increasing savings for large meshes. In general, reductions of up to 75% have been noted for some designs such as the Rutan Varieze. A test-suite of basic shapes and bodies has been tested to evaluate the fidelity of the new approach. This suite tested the current solver over thin wing plan-forms of various types. It was observed that the current theory is accurate for rectangular wings, 146 sweep wings, axisymmetric bodies and delta wings with all aerodynamic load predictions within 5-10% of experimental data at low angles of attack. At higher angles of attack, expected deviations between predictions and experiments are seen and these are attributed to the non-linearity of the physics close to, and after, flow separation. The validity of the Kutta condition is shown for the axisymmetric bodies when the vorticity shedding is eliminated and no trailing edges are marked. An advanced test-suite has also been developed to stress-test the solver as well as to test out the geometry handling pipeline required to handle these test cases. Through these tests, it was determined that the solver is able to provide high fidelity results for a wide variety of test cases. Results show a high degree of fidelity for typical fighter configuration geometries such as the F-16 and the F-18. In the case of both these geometries, the aerodynamic loads were found to be within ?10% of the experimentally obtained values. For the case of the F-16 it was shown that the range of applicability of the compressibility correction factor extends to Mach Numbers between 0.8-0.9. The F- 18 test case showed substantial improvement of numerical drag prediction with the current approach over contemporary solvers. Tests were also done for commercial airliner representative test cases. The DLR- F6 geometry was used and aerodynamic load predictions were within ?5-10% with slight over-predictions in drag and slight under-prediction in lift coefficient. This test was compared with CFD results using similar sized surface meshes as a starting point and the results showed substantially lower fidelity of solutions compared to experiments. In 147 addition, a reduction in CPU time was noticed from ~36 minutes (CFD) to about 3.5 minutes (current solver); a reduction of ~90%. Testing also included canard equipped geometries to determine the coupling effect of the lifting components in such geometries. Two designs were tested: the NASA Canard Fighter and the Rutan Varieze. In both cases the coupling effect of the canard on the main wing was observed under differing conditions and compared with experimental data. Results showed aerodynamic load predictions within ?5-10% and the coupling effect was readily visible through the use of the relaxed wake-strand model of vorticity shedding. Finally, the stated objective of this effort to develop a solver able to handle diverse geometrical cases hitherto restricted to Navier-Stokes solver has also been reasonably demonstrated by the various models and designs tested within the basic and advanced testing suites. Realistic designs such as the F-16, F-18 and Varieze geometries as well as reference cases such as the DLR-F6 (for CFD studies) and the NASA Canard Fighter (for legacy and contemporary potential-flow studies) illustrate the spread of designs that can be readily tested for aerodynamic performance using the Method of Integrated Circulation. REFERENCES 1. Reneaux, J., and Thibert, J-J., ?The Use of Numerical Optimization for Airfoil Design,? AIAA 3rd Applied Aerodynamics Conference, AIAA Paper 85-5026, Colorado Springs, CO, Oct. 1985 2. Ahuja, Vivek and Hartfield, Roy, ?Optimization of Scramjet Combustor designs using Genetic Algorithms and CFD?, AIAA 2008-5925 presented at the AIAA International Multidisciplinary Design and Optimization Conference (MDO), Victoria City, British Columbia, Canada, September 2008 3. Ahuja V., Hartfield R. J., ?Optimization of UAV Geometries for Aerodynamic Performance using Advanced Paneling Aero-Propulsive Paneling Models and Genetic Algorithms?, 6th AIAA Multidisciplinary Design Optimization Specialist Conference, Orlando, Florida , April 2010 4. Ahuja, V and Hartfield, Roy, ?Optimization of Air-breathing Hypersonic Aircraft Design using Euler Codes and Genetic Algorithms?, to be presented at the 46th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit , Nashville, Tennessee, USA, in July 2010 5. Sytsma, H.S., Hewitt, B.L. and Rubbert, P.E., ?A comparison of panel methods for subsonic flow computation,? AGARD AG-241, 1979 149 6. Wang, H.T., ?Comprehensive Evaluation of Six Thin-Wing Lifting Surface Computer Programs?, Naval Ship Research and Development Center, Rept. 4333, June 1974 7. Langan, T.J. and Wang, H. T., ?Evaluation of Lifting Surface Programs for Computing the Pressure Distribution on Planar Foils in Steady Motion?, Naval Ship Research and Development Center, Rept. 4021, May 1973 8. Maskew, B., ?Program VSAERO Theory Document?, NASA Contractor Report 4023, NASA-11945, September 1987 9. Maskew, B., ?Program VSAERO, A computer program for calculating the nonlinear aerodynamic characteristics of arbitrary configurations, User?s Manual,? NASA CR-166476, November 1982 10. Maskew, B., ?A surface vorticity method for calculating the pressure distribution over airfoils of arbitrary thickness and camber using two-dimensional potential flow,? Loughborough University of Technology report TT 6907, 1969 11. Epton, M.A. and Magnus A.E., ?PAN AIR- A computer program for predicting subsonic or supersonic linear potential flows about arbitrary configurations using a higher order panel method,? Version-I, Theory Document, NASA CR-3251, March, 1990 12. Ashby, D., Dudley, M., and Iguchi, S., ?Development and validation of an advanced low-order panel method,? NASA TM 101024, 1988 13. Dvorak, F.A., Woodward, F.A. and Maskew, B., ?A three-dimensional viscous/Potential flow interaction analysis methods for multi-element wings,? NASA CR-152012, July 1977 150 14. Katz, J. and Plotkin, A., ?Low-speed aerodynamics: from wing theory to panel methods,? McGraw-Hill Inc., ISB-0-07-050446-6, 1991 15. Maskew, B., ?Calculation of the three-dimensional potential flow around lifting non-planar wings and wing-bodies using a surface distribution of quadrilateral vortex rings,? Loughborough University of Technology, Dept. of Transport Tech. Report TT 7009, September 1970 16. Purcell, E.W., ?The vector method of solving simultaneous linear equations,? J. Math. Physics, Vol. 23, No. 180, 1953 17. Lamb, H., ?Hydrodynamics,? 6th ed., Cambridge University Press, 1932 18. Batchelor, G.K., ?Axial flow in trailing line vortices,? Journal of Fluid Mechanics, no. 20, pp. 645-658, 1964 19. Batchelor, G.K., ?An introduction to fluid mechanics,? Cambridge University Press, 1967 20. Filotas, L. T., ?Solution of the lifting line equation for twisted elliptical wings,? J. Aircraft, vol. 8, no. 10, pp. 835-936, 1971 21. Ahuja V., Hartfield R. J., ?Preliminary Drag Prediction Using Advanced Paneling Schemes?, 48th AIAA Aerospace Sciences Meeting, Orlando, Florida, January, 2010 22. Jobe, C.E., ?Prediction and Verification of Aerodynamic Drag, Part I:Prediction? Progress in Astronautics and Aeronautics, edited by C. E. Eugene, Vol. 98, AIAA, New York, 1985 23. Prandtl, L., ?Applications of modern hydrodynamics to aeronautics,? Report 116, NACA, 1921 151 24. Purvis, J.W., ?Simplified solution of the compressible lifting surface problem,? M.S. Thesis, Aerospace Engineering Department, Auburn University, Auburn, Alabama, August, 1975 25. Schlichting, H., and Truckenbrodt E., ?Aerodynamik des Flugzeuges,? vol. 2, chapter 7, Springer-Verlag, Berlin, 1969 26. Wentz, W.H. and Kohlman, D.L., ?Wind tunnel investigations of vortex breakdown on slender sharp-edged wings,? NASA Research Grant NGR-17-002- 043, University of Kansas Center for Research Inc., Report FRL 68-013, November 27, 1968 27. Polhamus, E.C., ?Application of the leading-edge suction analogy of vortex lift to the drag due to lift of sharp edge delta wings,? NASA TN D-4739, August 1968 28. Polhamus, E.C., ?A concept of the vortex lift of sharp-edge delta-wings based on a leading-edge suction analogy,? NASA TN D-3767, 1966 29. Wentz, W.H. and McMahon, M.C., ?Further experimental investigations of delta and double-delta wing flow fields at low speeds,? NASA CR-714, February, 1967 30. Weber J., Brebner G.G., ?Low speed tests on a 45-deg swept back wings, Part-I: pressure measurements on wings of aspect ratio 5?, ARC R&M 2882, 1958 31. Weissinger, J., ?The lift distribution of swept-back wings,? NACA TM-1120, 1947 32. Nguyen, Luat T. et.al, ?Simulator study of stall/post-stall characteristics of a fighter airplane with relaxed longitudinal static stability,? NASA TP-1538, December, 1979 152 33. Webb, T.S., Kent, D.R., ?Correlation of F-16 aerodynamics and performance predictions with early flight-test results,? AGARD conference proceedings, n 242, October 11-13, 1977 34. Banks, D.W., ?Wind tunnel investigation of the forebody aerodynamics of a vortex lift fighter configuration at high angles of attack,? SAW Paper 88-1419. October, 1988 35. Browne, L., Katz, J., ?Application of panel methods to wind-tunnel wall interference corrections,? 28th Aerospace Sciences meeting, Reno, Nevada, January 8-11, 1990 36. Levy, D. W., Zickuhr, T., Vassberg, J., Agrawal, S., Wahls, R. A., Pirzadeh, S., and Hemsch, M. J., ?Summary of Data from the First AIAA CFD Drag Prediction Workshop,? 40th AIAA Aerospace Sciences Meeting and Exhibit, AIAA Paper 2002-0841, Reno, NV, Jan. 2002 37. Laflin, K., Klausmeyer, S. M., Zickuhr, T.,Vassberg, J. C.,Wahls, R. A., Morrison, J. H., et al., ?Data Summary from Second AIAA Computational Fluid Dynamics Drag Prediction Workshop,? Journal of Aircraft, Vol. 42, No. 5, 2005, pp. 1165?1178 38. Vassberg, J. C., Tinoco, E. N., Mani, M., Brodersen, O. P., Eisfeld, B., Wahls, R. A., Morrison, J. H., Zickuhr, T., Laflin, K. R., and Mavriplis, D. J., ?Summary of the Third AIAA CFD Drag Prediction Workshop,? 45th AIAA Aerospace Sciences Meeting and Exhibit, AIAA Paper 2007-0260, Reno, NV, Jan. 2007 39. Stoll F., Koenig D.G., ?Large Scale Wind Tunnel Investigation of Close-Coupled Canard-Delta-Wing Fighter Model through high angles of attack?, AIAA Aircraft 153 Design, Systems and Technology Meeting, Fort Worth, Texas, 17-19 October 1983 40. Polhamus, E.C., Campbell, G.S., ?Aerodynamic characteristics of a wing with unswept quarter-chord line, aspect ratio 2, taper ratio 0.78 and NACA 65A004 airfoil section,? NACA Research Memorandum L50A18, March 08, 1950 41. Stoll F., Koenig D.G., ?Large Speed Wind Tunnel Measurements of a Canard Controlled Fighter Model through high angles of attack?, NASA TM-84403, 1983 42. Miller, S.G. and Youngblood, D. B., ?Application of USSAERO-3 and the PANAIR Production Code to the CDAF Model A Canard/Wing Configuration?, AIAA Paper 83-1829, 1983 43. Minter, E. A, Yates, R. W., ?Wind Tunnel Tests of a 0.4 Scale Fighter Model at High Angles of Attack-Analysis of Pressure Data?, NASA CR-166198, 1981 44. Levin, D., and Katz, J., ?A vortex-lattice method for the calculation of the non- steady separated flow over delta wings,? AIAA Paper 80-1803, 1980 45. Miller, S.G. and Youngblood, D.b., ?Application of USAERO-3 and the PANAIR Production code to the CADF Model A canard/wing configuration,?, AIAA paper 83-1829, 1983 46. Rutan, Burt, ?Development of a small high-aspect-ratio canard aircraft,? 1976 report to the Aerospace Profession, Tech Rev., vol. 13, no. 2, Soc. Exp. Test Pilots, 1976, pp-93-101 47. Long, P.Y., ?Wind tunnel investigation of a full-scale canard-configured general aviation airplane,? NASA TP-2382, March, 1985 154 48. Karamcheti, K., ?Principles of Ideal-Fluid Aerodynamics,? R.E. Krieger, Malabar, FL, 1980 155 APPENDIX-A DERIVATION OF THE BIO-SAVART LAW FOR VORTEX RINGS The continuity equation for incompressible and inviscid equation can be written as: ? = 0 (A.1) The velocity field can be written as the curl of a vector field B, such that: ? = ? (A.2) Since the curl of a gradient vector is zero, B is indeterminate to within the gradient of a scalar function of position and time, and the vector field B can selected such that: ? = 0 (A.3) The vorticity then becomes: ? = ? = ( ?) = ( ?) ? (A.4) Substituting Equation-A.3 into Equation-A.4, ? = ? (A.5) The solution to Equation-A.5 is arrived at using Green?s theorem by Karamcheti 48 and is written as: 156 ? = 14 ? ?| |? (A.6) Combining Equation-A.2 with Equation-A.6 gives: ? = 14 ? ?| |? (A.7) Since the infinitesimal volume can be split into an infinitesimal area normal to the Vorticity and an incremental distance ?? on the filament so that: ? = ? ?? (A.8) Also, circulation on the infinitesimal area ? is given as: = ?? (A.9) Combining Equation-A.8, Equation-A.9 and Equation-A.7, we get: ? = 14 ? | |?? = 4 ??? | | (A.10) Equation-A.10 is the required vector formulation of the Biot-Savart law as used in the current effort.