> ,'t " "" --" NASA Contractor Report 189739 H=P ADAPTIVE METHODS FOR FINITE ELEMENT ANALYSIS OF AEROTHERMAL LOADS IN HIGH-SPEED FLOWS H.J. Chang, J.M. Bass, W. Tworzydlo, and J.T. Oden COMPUTATIONAL MECHANICS COMPANY, INC. Austin, Texas Contract NAS 1-18746 January 1993 (N ASA- OR- I'-9_7'.3-7) H-P ADAPTIVE N93-18093 t_,ETHLIgS FOR FI=NITE ELEMENT ANALYSIS CP AC_i_-_THERMAL LOADS IN HIGH-SPEED r-L .j_ 3 (Computational Mechanics Uncl as Co. ) 290 p G3/34 0142862 National Aeronautics and Space Administration Langley Research Center Hampton, Virginia 23665-5225 Contents 1 Introduction 4 2 Formulation of the Governing Equations and Numerical Algorithms 6 2.1 Notation and Formulation of the Equations ................... 6 2.1.1 Nondimensional Form of the Navier-Stokes Equations ......... 8 2.2 Taylor-Galerkin Algorithms ........................... 10 2.2.1 One-Step Taylor-Galerkin Algorithm .................. 11 2.2.2 The Two-Step Algorithm ......................... 14 2.3 Boundary Conditions ............................... 18 2.3.1 Boundary Condition for the One--Step Algorithm ........... 19 2.3.2 Boundary Conditions for the Two-Step Method ............ 28 2.4 Artificial Viscosity ................................ 31 3 An h-p Finite Element Method 34 3.1 A Variational Formulation ............................ 35 3.2 Finite Element Approximation .......................... 36 3.3 Adaptivity ..................................... 37 3.4 Regular and Irregular Meshes 37 3.5 Basic Assumptions ................................ 40 3.6 Definition of an Element ............................. 40 3.7 Continuity for Regular Meshes .......................... 46 3.8 Continuity for 1-Irregular Meshes. Constrained Approximation ........ 48 3.9 Calculation of the Element Load Vector and Stiffness Matrix ......... 57 3.10 Constraints in the One-Dimensional Case ................... 60 3.11 Constraints for Two-Dimensional Subparametric Elements .......... 64 3.12 Constrained Approximation in a Three-Dimensional Case ........... 64 3.13 Constraints for a Wall .............................. 67 3.14 Constraints for an Edge ............................. 70 3.15 The Constrained Base Shape Functions ..................... 72 3.16 Interpretation of _;. Calculation of the Load Vector and Stiffness Matrix . . 72 111 PRECEDING PAGE B_|ANI!_ NOT FILMED 3.17 Concluding Remarks on Constrained Approximation ............. 74 3.18 Some Details Concerning the Data Structure .................. 74 4 Adaptivity 79 4.1 Error Estimation Techniques ........................... 79 4.1.1 Interpolation Error Estimate ....................... 79 4.1.2 Residual Error Estimate ......................... 80 4.1.3 Relative Error Estimate ......................... 93 4.2 Directional Adaptation Indicator ........................ 102 4.3 Adaptive Strategies ................................ 114 5 Implicit/Explicit Procedures 115 5.1 Formulation of implicit/explicit schemes .................... 115 5.2 Selection of implicit and explicit zones ..................... 117 5.3 Computational procedure ............................ 121 6 Some Nonstandard Algorithms 133 6.1 Integration and Underintegration Procedures .................. 133 6.2 Routines for Performing Refinements in Boundary Layer Zone ........ 137 6.3 Postprocessing in Three Dimensions ....................... 137 6.4 Solution Correcting Procedures ......................... 139 7 Numerical Examples 140 8 Phase II Project Summary and Future Directions 268 9 References 270 A Performance Issues 273 A.1 Iterative Methods for the h-p Finite Element Method ............. 273 A.I.I Basic Description of the Iterative Method ............... 274 A.I.2 Separating a Preconditioner From an Accelerator ........... 275 A.I.3 Patch Definition ................... . .......... 276 A.I.4 Numerical Results ............................ 278 iv A.2 Matrix Condensation ............................... 283 V Project Summary Over the past four years the Computational Mechanics Company, Inc. has been the principal investigator into the development of a new class of computational methods for modeling hypersonic viscous flows in two- and three-dimensional domains. The ultimate goal of this project is to provide NASA with a highly accurate computational tool for resolving fine scale flow features typical of shock wave interactions and strong viscous interaction regions in hypersonic flows. Toward this end, a research and development effort was put forth in the area of adaptive computational finite element methods for high speed flows based on unstructured mesh concepts and employing local estimates of the solution error to optimally change the computational grid to minimize the numerical error. The approach developed here, which combines h-adaptive and anisotropic p-enrichment mesh modification procedures, has proven for certain classes of problems to provide exponential rates of convergence of the solution using possibly an order of magnitude fewer degrees of freedom than conventional methods while obtaining the same level of computational accuracy. During the first phase of this project (years 1-3), work was focused on a number of research topics which were crucial to the success of h-p adaptive finite element methods for hypersonic flows. Many of these issues were resolved and the results are presented again in detail in the body of this report for completeness. Summarizing some of the more significant developments of this first phase of the effort: • The development of the first h-p finite element data structure for quadrilateral and hexahedral elements. • The development of an h-p adaptive package that includes m A local error estimation capability for driving the adaptive strategy. A spectral enrichment and h-refinement methodology to change the structure of the computational grid. The formulation of a generalized methodology for handling nodal point constraints for higher order polynomials. (Note that this difficulty does not occur with either h-refinement or p-enrichment individually.) • The development of an algorithm for manual anisotropic enrichment of both two- dimensional and three-dimensional elements. • The formulation and implementation of a version of a preconditioned block Jacobi- GMRES method for higher order spectral elements. • The formulation andimplementationofa one-stepTaylor-Galerkinsolution algorithm for h-p methods. The formulation and implementation of an implicit two-step Taylor-Galerkin solution algorithm which solves first an Euler step followed by a viscous step. An investigation of artificial dissipation mechanisms appropriate for h-p adaptive com- putational methods with possibly highly distorted elements or elements with high as- pect ratios. In addition to these efforts on research-oriented topics, a user friendly, graphics oriented, interface for the two- and three-dimensional codes was developed. This interface includes both a batch and interactive option and full graphics capabilities for displaying the solu- tions, extracted quantities, and plots of the local estimates of the computational error. A simple grid file interface was also developed to read in neutral grid files generated either by the GAMMA2D or GAMMA3D codes, developed in-house at COMCO, or by PATRAN. (The formats necessary for PATRAN and the neutral file interface are discussed in the grid generation section of the user manual.) The second phase of this effort, conducted over the last 12 months, has focused on two special research and development topics which are in general related to the performance of the flow solver. These topics include: • The development of implicit/explicit computational methods (in two dimensions) for integrating the Navier Stokes equations forward in time. • The development of computational methodologies and algorithms which provide for automated directional p-enrichment of the computational mesh. A third topic on which we dedicated considerable computer resources was a continued investigation of artificial dissipation mechanism for h-p adaptive computational methods. In particular, numerous test cases were run for the Mach 14 Holden problem using various artificial mechanism to determine an appropriate model for capturing the recirculation bubble and other fine features of the solution. All of the capabilities and options developed during the first and second phases of the effort have been implemented and/or tested in a two-dimensional and/or a three-dimensional finite element code. Using these codes a number of test cases have been solved to verify the functionality of the software. These benchmark cases include flow past a blunt body with an incident shock to produce a Type IV (Edney) interaction, a Holden problem with inflow Mach number 14, and a rearward facing step with a strong expansion. Experimental data 2 on each of these benchmark problems is available and comparisons with this data are made throughout the results section of this report. 1 Introduction The commitment to develop the National Aero-Space Plane and Maneuvering Reentry Ve- hicles has generated resurgent interest in the technology required to design structures for hypersonic flight. As these vehicles cruise, accelerate and/or decelerate in the atmosphere, highly complex patterns of shock wave interactions and shock wave boundary layer inter- actions develop which produce severe local pressures and extreme local heating rates. To provide adequate safety factors in the thermal-structural designs an accurate determination of the aerothermal loads, especially in the local areas of strong shock interactions and strong viscous interactions such as interacting shock waves on a leading edge, shock/boundary layer interactions, corner flows, etc., is required. In general, ground based test facilities can pro- vide only limited data for the expected flight conditions, at a considerable cost, and thus designers must depend heavily on analytically predicted aerothermal loads. The analyti- cal methodologies for these predictions must be accurate (to within some predetermined measure), robust, computationally economical, and geometry independent. For a designer, the prediction of aerothermal loads of hypersonic structures is one of the most challenging problems in computational fluid dynamics. It requires flow analysis and shock prediction at very high Mach numbers, as well as a realistic calculation of aerother- mal loads. These tasks required, in general, advanced computational strategies, extremely accurate discretization techniques, and powerful postprocessing capabilities. During the first two phases of this project we have developed, at the Computational Mechanics Company, an implicit/explicit, anisotropic h-p finite element methodology for the analysis of high speed flows and the prediction of aerothermal loads. The basic idea of the h-p version of the finite element method is to combine local mesh refinement (an h-method) with anisotropic polynomial enrichment (a p-method) in order to achieve con- vergence rates not attainable with fixed mesh methods or with any of the above methods applied separately. In practical terms, the h-p method means maximum numerical accuracy at a minimal computational cost. The remainder of this report summarizes the results of the first and second phases of a four year research and development effort oriented toward the resolution of several theoretical and computational issues related to the above problems. In Section 2 the basic formulation of the Navier-Stokes equations which governs the compressible viscous flow is presented. This is followed by a detailed discussion of a general family of implicit Taylor Galerkin algo- rithms for solving the compressible flow equations, and a review of the numerical boundary conditions associated with the algorithms. The final part of this section provides a brief discussion about artificial dissipation mechanisms which may be more appropriate for highly distorted elements or elements with high aspect ratios. The Taylor-Galerkin formulation 4 and the general boundary condition treatment presented in Section 2 is excerpted from our published work on these results (e.g. [8, 10, 31]). The next section, Section 3, provides a detailed review of the h-p finite element methodology. Included in this review is a discussion of the data structure, hierarchical shape functions, constrained approximation, and assump- tions and restrictions which have been incorporated into the h-p formulation. The work reported in Sections 3.3-3.11 is based on extensions of the general h-p adaptivity methods of Oden, Demkowicz, and Rachowicz, originally published in 1989 (see [9, 25]) and expanded and generalized to apply to thermomechanical load predictions in this project. Section 4 follows with a discussion of the adaptive strategy which is used in conjunction with the mesh modification algorithms. Here the error estimates used to drive the adaptive package are outlined and the procedure for adapting the computational mesh is presented. Section 5 follows with a detailed discussion of the implicit/explicit methodologies that have been implemented within the context of the two-dimensional code. In particular, in this section we review the basic concepts associated with using implicit/explicit methods and how one would select implicit end explicits zones within the computational domain. This section con- cludes with a brief overview of other computational procedures also required to efficiently implement an implicit/explicit methodology. The contents of this report includes a summary of the efforts completed during both the first and second phases of the development effort. During the second phase of the project, our efforts have focused on implicit/explicit methods, directional enrichment techniques and algorithms, and optimal artificial dissipation mechanisms. The details of this effort are provided in sections 2.2, 4.2, 5, and 7. The other basic sections are essentially the same as in previous reports except for various corrections. The next section, Section 6, presents some non-standard algorithms used in the element calculations and the postprocessing module. This is followed by Section 7, which presents the results of the numerical test cases that have been run over the course of this project. This includes both simple test problems used for verification purposes as well as the benchmark problems supplied by NASA-Langley. The final section provides a brief discussion of some possible directions for future research and development in the area of h-p adaptive methods for hypersonic flows. 5