Citation
On total variation diminishing schemes for transonic turbulent flow computation

Material Information

Title:
On total variation diminishing schemes for transonic turbulent flow computation
Creator:
Chen, Ming-Hsiung, 1951-
Publication Date:
Language:
English
Physical Description:
xi, 138 leaves : ill. ; 28 cm.

Subjects

Subjects / Keywords:
Approximation ( jstor )
Conservation laws ( jstor )
Entropy ( jstor )
Jacobians ( jstor )
Matrices ( jstor )
Pressure ( jstor )
Projectiles ( jstor )
Scalars ( jstor )
Trucks ( jstor )
TVD schemes ( jstor )
Aerodynamics, Transonic ( lcsh )
Aerospace Engineering, Mechanics, and Engineering Science thesis Ph. D
Computational fluid dynamics ( lcsh )
Dissertations, Academic -- Aerospace Engineering, Mechanics, and Engineering Science -- UF
Navier-Stokes equations ( lcsh )
Turbulence ( lcsh )
Genre:
bibliography ( marcgt )
non-fiction ( marcgt )

Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 1988.
Bibliography:
Includes bibliographical references.
General Note:
Typescript.
General Note:
Vita.
Statement of Responsibility:
by Ming-Hsiung Chen.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
Copyright [name of dissertation author]. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Resource Identifier:
024931662 ( ALEPH )
20117342 ( OCLC )
AFM5797 ( NOTIS )
AA00004795_00001 ( sobekcm )

Downloads

This item has the following downloads:


Full Text
xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID EOKPRX157_XJMVF0 INGEST_TIME 2011-11-21T19:37:35Z PACKAGE AA00004795_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES



PAGE 1

ON TOTAL VARIATION DIMINISHING SCHEMES FOR TRANSONIC TURBULENT FLOW COMPUTATION BY MING-HSIUNG i CHEN A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1988 ILPJJ LIBRARIES

PAGE 2

ACKNOWLEDGEMENTS The author wishes to express his sincere appreciation to his thesis advisor, Dr. Chen-Chi Hsu, for his guidance and understanding throughout this endeavor. Special appreciation is expressed to Dr. W. Shyy for his spending countless hours reading the manuscripts and suggesting improvements in this thesis. Also special thanks go to other thesis committee members, Dr. U.H. Kurzweg, Dr. L.E. Malvern, and Dr. A.K. Varma for their time and comments. The author wishes to thank Dr. H.C. Yee of NASA Ames Research Center for providing a listing of her 2-D TVD code, which was very helpful for the development of the current research code for axisymmetric aerodynamic problems . All the numerical calculations reported in this thesis were performed on CRAY XMP/48 supercomputer at the National Center for Supercomputing Applications (NCSA) , University of Illinois at Urbana-Champaign. Finally, the deepest gratitude is expressed to his wife Ko-Hui for her support and encouragement. ii

PAGE 3

TABLE OF CONTENTS PAGE ACKNOWLEDGEMENTS ii LIST OF TABLES V LIST OF FIGURES vi ABSTRACT X I . INTRODUCTION 1 II . PRELIMINARY 12 2 . 1 Hyperbolic Conservation Laws 12 2 . 2 Background of TVD Schemes 20 III . UPWIND AND SYMMETRIC TVD SCHEMES 25 3 . 1 First Order Accurate TVD Scheme 25 3 . 2 Second Order Accurate TVD Scheme 30 3.3 Extension of Scalar TVD Schemes to a System of Conservation Laws 33 3 . 4 Upwind TVD Schemes 39 3 . 5 Symmetric TVD Schemes 40 3.6 Linearized Conservative Implicit (LCI) Form... 42 3.7 Alternative Direction Implicit (ADI) Form 44 IV. FLUX LIMITERS 46 V. ON THREE-DIMENSIONAL SYMMETRIC AND UPWIND TVD SCHEMES IN GENERALIZED COORDINATES 63 5.1 The 3-D Equation in Generalized Coordinates. .. 63 5.2 Algorithm of 3-D Schemes in Generalized Coordinates 72 5.3 A Conservative Linearized ADI Form for Steady-State Appliciation 75 VI . THIN LAYER NAVIER-STOKES SOLVERS 78 6 . 1 Thin-Layer Approximation 79 6.2 Azimuthal (or r?-) Invariant Equations 80 6 . 3 Turbulent Model 81 6.4 Implicit Central Difference Algorithm !.'.'s2 iii

PAGE 4

6.5 Dissipation Model of TVD scheme 84 6.6 Grid System 85 6.7 Numerical Boundary Conditions 88 6 . 8 Spatial-Varying Time Step 89 VII. RESULTS AND DISCUSSION 91 7 . 1 Numerical Implementation 91 7.2 Result of Ma, = 0.96 94 7.3 Grid Density and Numerical Accuracy 101 7 . 4 Result of M*, = 0.91 112 7 . 5 Result of M*, = 1.2 117 VIII . SUMMARY AND CONCLUDING REMARKS 122 APPENDIX A: NAVIER-STOKES EQUATIONS AND THIN LAYER APPROXIMATION 124 APPENDIX B: BEAM-WARMING BLOCK ADI ALGORITHM AND ARTIFICIAL DISSIPATION 130 REFERENCES 1 34 BIOGRAPHICAL SKETCH 138 IV

PAGE 5

LIST OF TABLES PAGE Table 7.1 Pressure and Shear Drags for Different Scheme of 90x60 Grid at M«, = 0.96 loo Table 7.2a L 2 -Norm Residuals and CPU Time for Different Schemes of 90x60 Grid at M^ = 0.96 Using Fixed Time Step 101 Table 7.2b L 2 -Norm Residuals and CPU Time for Different Schemes of 90x60 Grid at M« = 0.96 Using Fixed Time Step 101 Table 7.3 L 2 -Norm Residuals and CPU Time for Different Schemes of 70x40 Grid at M*. = 0.96 112 Table 7.4 Pressure and Shear Drag for Different Schemes of 70x40 Grid at K^ = 0.96 112 Table 7.5 L 2 -Norm Residuals and CPU Time for Different Schemes of 70x40 Grid at M,, = 0.91 115 Table 7.6 Pressure and Shear Drag for Different Schemes of 70x40 Grid at M^ =0.91 115 Table 7.7 Pressure and Shear Drag for Different Schemes of 70x40 Grid at M^ = 1.2 121

PAGE 6

Figure 4 . 1 Figure 4.2 Figure 4 . 3 Figure 4 . 4 Figure 4 . 5 Figure 4 . 6 Figure 6 . 1 Figure 6.2 Figure 7 . 1 Figure 7.2a Figure 7.2b Figure 7.2c Figure 7 . 2d LIST OF FIGURES PAGE TVD region 51 Second order TVD region 52 Van Leer's limiter, Eq.(4.17) 54 Minmod limiter, Eq. (4 . 18) 56 Roe's superbee limiter, Eq.(4.19) 58 Graphical representation of three limiters, Eq. (4.6a,b,c) 62 The configuration of SOCBT projectile 86 Mapping from physical domain to computational domain 87 The 90x60 "C-type" grid for the SOCBT projectile with sting 93 The surface pressure coefficients of Beam-Warming scheme for the SOCBT projectile M*, = 0.96, and Re = 7.7xl0 5 93 The surface pressure coefficients of upwind TVD scheme with Limitl for the projectile with M*, = 0.96, and Re = 7.7xl0 5 using a 90x60 C grid 96 The surface pressure coefficients of upwind TVD scheme with Limit2 for the projectile with M*, = 0.96, and Re = 7.7xl0 5 using a 90x60 C grid 96 The surface pressure coefficients of symmetric TVD scheme with Limit3 for the projectile with M,, = 0.96, and Re = 7.7xl0 5 using a 90x60 C grid 97 vi

PAGE 7

Figure 7.3 Figure 7.4a Figure 7.4b Figure 7.4c Figure 7.4d Figure 7.5a Figure 7.5b Convergence history of different schemes (a) Beam-Warming scheme, (b) Upwind TVD scheme with Limit2, (c) Upwind TVD scheme with Limitl and (d) Symmetric TVD scheme with Limit3 The surface pressure coefficients of Beam-Warming scheme for projectile with M« = 0.96, and Re = 7.7xl0 5 using a 60x60 C grid. The L 2 -norm residual of 5.16xl0~ 5 can be reached in around 5000 time steps 99 103 The surface pressure coefficients of upwind TVD scheme with Limitl for projectile with M« = 0.96, and Re = 7.7xl0 5 using a 60x60 C grid. The L2~norm residual of 3.9xl0~ 5 can be reached in around 1500 time Steps 103 The surface pressure coefficients of upwind TVD scheme with Limit2 for projectile with M« = 0.96, and Re = 7.7xl0 5 using a 60x60 C grid. The L2~norm residual of 8.5xl0~ 5 can be reached in around 1900 time steps 104 The surface pressure coefficients of symmetric TVD scheme with Limit3 for projectile with M« = 0.96, and Re = 7.7xl0 5 using a 60x60 C grid. The L 2 -norm residual of 3.4xl0" 5 can be reached in around 1100 time steps 104 The surface pressure coefficients of Beam-Warming scheme for projectile with M« = 0.96, and Re = 7.7xl0 5 using a 90x40 C grid. The L 2 -norm residual of 2.59xl0" 6 can be reached in around 3000 time steps 105 The surface pressure coefficients of symmetric TVD scheme with Limit3 for projectile with M«. = 0.96. and Re = 7.7xl0 5 using a 90x40 C grid. The L 2 -norm residual of 4.8xl0~ 5 can be reached in around 1100 time steps 105 Vll

PAGE 8

Figure 7 . 6a Figure 7.6b Figure 7.6c Figure 7.7a Figure 7.7b Figure 7.7c Figure 7.8 Figure 7.9a The surface pressure coefficients of upwind TVD scheme with Limitl for projectile with M« = 0.96, and Re = 7.7xl0 5 using an 80x40 C grid. The L 2 -norm residual of 7.68xl0" 6 can be reached in around 1300 time steps.. 107 The surface pressure coefficients of upwind TVD scheme with Limit2 for projectile with M» = 0.96, and Re = 7.7X10 5 using an 80x40 C grid. The L 2 -norm residual of 8.15xl0~ 6 can be reached in around 1500 time steps.. 107 The surface pressure coefficients of symmetric TVD scheme with Limit3 for projectile with M« = 0.96, and Re = 7.7xl0 5 using an 80x40 C grid. The L 2 -norm residual of 5.4xl0~ 5 can be reached in around 1100 time steps. The surface pressure coefficient of upwind TVD scheme with Limitl for projectile with M*. = 0.96, and Re = 7.7xl0 5 using a 70x40 C grid 108 109 The surface pressure coefficient of upwind TVD scheme with Limit2 for Projectile with H*, = 0.96, and Re = 7.7xl0 5 using a 70x40 C grid... 109 The surface pressure coefficient of symmetric TVD scheme with Limit3 for projectile with M» = 0.96, and Re = 7.7xl0 5 using a 70x40 C grid no The convergence history of different TVD schemes at M^ =0.96 using a 70x40 grid, (a) Upwind TVD scheme with Limit2, (b) Upwind TVD scheme with Limitl and (c) Symmetric TVD scheme with Limit3 Ill The surface pressure coefficients of upwind TVD scheme with Limitl for projectile with M*, = 0.91, and Re = 7.7xl0 5 using a 70x40 C grid 113 Figure 7.9b The surface pressure coefficient of upwind TVD scheme with Limit2 for projectile with M^ = 0.91, and Re = 7.7xl0 5 using a 70x40 C grid... 113 Vlll

PAGE 9

Figure 7 . 9c Figure 7.10 Figure 7.11a Figure 7.11b Figure 7.11c Figure 7 . 12 The surface pressure coefficients of symmetric TVD scheme with Limit3 for projectile with M*, = 0.91, and Re = 7.7xl0 5 using a 70x40 C grid 114 The convergence history of different TVD schemes at M^ = 0.91 using a 70x40 grid, (a) upwind TVD scheme with Limit2, (b) upwind TVD scheme with Limitl and (c) symmetric TVD scheme with Limit3 116 The surface pressure coefficients of upwind TVD scheme with Limitl for projectile with M,,, = 1.2, and Re = 7.7xl0 5 using a 70x40 C grid.... 118 The surface pressure coefficients of upwind TVD scheme with Limit2 for projectile with M^ = 1.2, and Re = 7.7xl0 5 using a 70x40 C grid 118 The surface pressure coefficients of symmetric TVD scheme with Limit3 for projectile with M*, = 1.2, and Re = 7.7x10 s using a 70x40 C grid 119 The convergence history of different TVD schemes at H, = 1.2 using a 70x40 grid, (a) upwind TVD scheme with Limit 2, (b) upwind TVD scheme with Limitl and (c) symmetric TVD scheme with Limit3 120 IX

PAGE 10

Abstract of Dissertation Presented to Graduate School of the University of Florida in Partial Fulfillment of the Requirement for the Degree of Doctor of Philosophy ON TOTAL VARIATION DIMINISHING SCHEMES FOR TRANSONIC TURBULENT FLOW COMPUTATION By MING-HSIUNG CHEN December, 1988 Chairman : Chen-Chi Hsu Major Department : Aerospace Engineering, Mechanics & Engineering Science Recently a number of new techniques for constructing nonlinear, second-order accurate, high-resolution shockcapturing schemes for systems of hyperbolic conservation laws have been developed in order to simulate complex flow fields more accurately and efficiently. These schemes, called total variation diminishing (TVD) , are free from generating spurious oscillations across the shocks and contact discontinuities, and can converge to a physically realizable solution with the aid of entropy correction. The extension work has been conducted to extend the TVD schemes developed for two-dimensional planar flows to

PAGE 11

three-dimensional flows. On the application side, a thinlayer Navier-Stokes flow code based on the conventional Beam-Warming scheme has been modified to incorporate the TVD schemes in the solution process. A variable time step procedure has also been implemented into the code for more efficient computation of the steady state flows. The transonic turbulent flow over an axisymmetric projectile is adopted as the test problem to assess the performance of three variants of the TVD scheme, two upwind schemes with flux limiters proposed by Van Leer and Roe, and one symmetric scheme. The results obtained show that all three variants can accurately resolve the shock and flow profiles with fewer grid points than the Beam-Warming scheme. Furthermore, although the computing cost of the TVD schemes is higher than the Beam-Warming scheme for each time step, the higher convergence rate of these schemes makes them computationally more efficient overall. The combinations of the high accuracy, good robustness, and improved computational efficiency offered by the TVD schemes make them extremely attractive for computing flow with shocks. Among the three variants, it is found that the convergence rate of the flux limiter proposed by Van Leer for the upwind TVD scheme is least affected by the incoming flow Mach numbers. xi

PAGE 12

CHAPTER I INTRODUCTION Computational fluid dynamics (CFD) has recently become an increasingly powerful tool in the aerodynamic design of aerospace systems. A number of CFD codes have been developed to simulate the high-speed flow problem. Routine calculations can now be performed to partially substitute for the wind tunnel measurements. Consequently, designs can now be analyzed and improved in shorter time and with lower cost. However, despite these fast advancements, much remains to be done to improve both the accuracy and efficiency of the numerical flow calculations. One major topic particularly relevant to aerospace applications is the treatment of shocks. In general there are three major difficulties in computing flow with shocks. These difficulties are listed as follows: (i) Schemes that are second (or higher) order accurate, such as the well-known methods developed by Lax and Wendroff [1], MacCormack [2], and Beam and Warming [3], may produce spurious oscillations across the discontinuities. First order schemes, on the other hand,

PAGE 13

2 suffer from severe numerical diffusion. They cannot provide accurate solutions for complicated flow fields unless an excessively fine grid distribution is adopted. (ii) The numerical schemes may develop nonlinear instabilities when discontinuities are encountered. (iii) The schemes may yield nonphysical solutions. It is well known that these classical schemes are not robust and accurate enough for flow containing discontinuities . In order to suppress high-freguency numerical errors of second or high order schemes associated with the shock calculation, the most common procedure is to add artificial dissipation terms. For example, in the widely used BeamWarming scheme, a fourth-order explicit dissipation is used to control non-linear instabilities, whereas the secondorder implicit dissipation is built-in to relax the stability bound from the explicit fourth-order dissipation. All these numerical dissipation terms contain adjustable parameters. For inviscid steady-state airfoil calculation, the use of the above type of dissipation terms may still produce spurious oscillations near shocks [4]. Such oscillations are not only undesirable but often lead to a breakdown in stability of the numerical scheme. This is the major drawback of the Beam-Warming scheme and other schemes .

PAGE 14

3 In order to overcome these difficulties, a number of new techniques for construction of nonlinear, second-order accurate, high-resolution schemes for hyperbolic conservation laws have been developed in recent years in order to simulate complex flow fields more accurately and efficiently [5]. All these schemes share the following attractive properties: (1) The schemes do not generate spurious oscillations. They are total variation diminishing in the nonlinear scalar case and the constant coefficient system case [5]. (2) They are consistent with the conservation law form and entropy inequality [6]. This property guarantees that the numerical weak solutions of the hyperbolic conservation law are physically realizable. In the following, a brief review is conducted to trace the original ideas that the modern TVD schemes are based on, and to present several variants of the TVD scheme proposed in the literature which have proven to be useful. CIR Method The earliest method for gas-dynamic equations in characteristic form was proposed by Courant et al. [7]. Their procedure, sometimes called the CIR (CourantIssacson-Rees) method is to trace back from (iAx, (n+l)At) all three characteristic paths. Since the problem is nonlinear, the directions of these paths are not known

PAGE 15

exactly, but to a first approximation they can be taken equal to their known directions at (iAx, nAt) . Then each characteristic equation is solved by using interpolated data at time nAt, in the interval to the left of i for characteristics with positive speed, and in the interval to the right of i for characteristics with negative speed. This method has two principal drawback: because it does not deal with the integral form (conservative form) of the equations, it cannot convey a shock wave with the proper speed; also, the resulting method is only first-order accurate. Godunov's scheme In 1959, Godunov [8] proposed the idea of advancing the solution to the next time-level by solving a set of Riemann problems. It is noted that the Riemann problem for any hyperbolic system of conservation laws arises if initial data are prescribed as two semi-infinite states (u = u L for x < 0, u= u R for x > ) For one-dimensional Euler equations, the solutions consist of three waves (u, u±c) which represent the propagation speed of the contact discontinuity, expansion wave, and shock respectively. To describe Godunov's use of the Riemann problem, let u j+l/2 °e the average state over (j±l/2)Ax at (nAt). He then replaced the data by an approximate distribution in which the state inside each interval is uniform and equal

PAGE 16

5 to u j . For each interface (j+l/2)Ax one can solve the Riemann problem with ur = uj+^ and ul = u j . This gives an exact solution Uj+1/2 to tne approximate problem, assuming that At is small enough that the waves do not interact from neighboring interfaces. The solution at (n+l)At can again be approximated by a piecewise uniform distribution, and then the process can be repeated. The advantage of Godunov's scheme is the clear physical picture of interaction involved, and the absence of any arbitrary parameter. Even though Godunov's scheme is monotonicity-preserving and produces good steady shock profile, the main drawback of Godunov's scheme are that it is cumbersome to compute and that is involves the high cost to solve the Riemann problem exactly. Also this scheme is only first -order accurate. Flux Correct Transport (FCT) Scheme The flux Correct Transport (FCT) scheme was first proposed by Boris and Books in 1973 [9], and later extended and refined by Zalesak [10]. The FCT scheme may be considered as a two-step method, a lower order transport stage corrected by a higher order antidif fusive flux in such a way as to prevent the occurrence of spurious oscillations. In general, The FCT scheme performs well compared to other schemes. However, as Zalesak has pointed out, peaked data clipping can occur, and he proposed a

PAGE 17

6 modified flux limiter which eliminates the clipping and is easy to extend to multidimensions. Harten's Scheme In 1983, Harten [5] developed a simple algebraic criterion for constructing a first-order accurate TVD scheme. He then added a limited anti-diffusive flux to a first order accurate entropy-satisfying scheme to obtain the second-order accurate TVD scheme. This technique is referred as the 'modified flux approach'. The modified flux is chosen so that the scheme is second-order at a smooth region and will switch itself into first-order accuracy at the vicinity of points of extrema. In [11] he also has extended the class of explicit TVD schemes to a more general category which includes a one-parameter family of implicit second-order TVD schemes. More recently Harten's TVD scheme was modified and generalized by Yee [12] and implemented to solve the two-dimensional Euler equations of gas dynamics for the airfoil problem. The numerical results indicate that Harten's scheme is both robust and accurate. MUSCL ( Monotonic Upwind Schemes of Conservation Laws ) Schemes B. van Leer [13] adopted an approach developed by Godunov (1959), who used a control volume formation to

PAGE 18

7 ensure that the scheme would be conservative and therefore dealt with averaged values instead of nodal values, Van Leer also observed that it is beneficial to work with averaged gradients in addition to averaged values, i.e., one can obtain second-order accuracy in Godunov's scheme by replacing the piecewise-constant initial data of the Riemann Problem with piecewise-linear initial data. The slope of the piecewise-linear initial data is selected so that spurious oscillation will be prevented. His result indicated that one does not have to resort to a Godunov type of Riemann solver in order to obtain good shock resolution. Colella and Woodward [14] introduced the piecewise parabolic method (PPM) ; They further refined van Leer's idea by using piecewise parabolic initial data, and the results they obtained are remarkably good. Roe's Scheme Roe [15,16,17] has developed a numerical fluctuation approach toward schemes for the numerical solution of both scalar hyperbolic conservation law and hyperbolic systems of conservation laws. His approach combines a novel interpretation of difference schemes with higher order accuracy and monotonicity preservation and also incorporates these ideas into an approximate Riemann solver for hyperbolic systems. Roe's approximate Riemann solver requires that there exists an average function A(u R ,u L )

PAGE 19

8 which is a nonlinear function of up and ul. The average function is constructed by the use of flux limiters. Using this approach, Roe has obtained solutions with a crisp shock, free from the spurious oscillations inherent to classical shock-capturing schemes. Unfortunately, Roe's scheme may admit nonphysical (entropy violating) solutions, i.e., expansion shocks. Various remedies have been proposed [18], [19]/ [20]. Osher's Scheme Osher's scheme [21] is based on an approximate Riemann solver, but compression and rarefaction waves are used to approximate shocks. This leads to a simpler and smoother algorithm. The numerical flux functions, which are at least continuously differentiable, are written in closed form and include various switches which make them upwind. Also the limit solutions satisfy the entropy condition. It is observed that Osher's scheme requires more operations than Harten's or van Leer's [22]. Symmetric TVD Scheme In 1984, Davis [23] showed that a certain class of TVD schemes can be interpreted as a Lax-Wendroff scheme plus and upwind weighted conservative numerical dissipation term. He then simplified the scheme by eliminating the upwind weight of this numerical dissipation term and also

PAGE 20

9 ensured that the simplified scheme still possessed the TVD properties. Shortly after that, Roe [24] reformulated Davis's scheme in such a way that the resulting scheme was easier to analyze and included a class of TVD schemes not obtained by Davis. More recently Yee [25] generalized Roe • s schemes to a one-parameter family of second-order explicit and implicit TVD schemes. Therefore, the formulation of Roe and Davis ' s scheme can simply be viewed as special cases of Yee's explicit symmetric TVD schemes. Also the main advantages of Yee's formulation are that stiff problems can be solved more efficiently by using implicit method and that steady-state solutions are independent of the time-step. Present Contributions In the present study, the viewpoint of interpreting the TVD scheme as a modified flux is adopted in order to make the task of practical implementation more tractable. One can modify a standard three-point central difference code ( e.g. Beam-Warming algorithm ) by simply changing the conventional dissipation terms into the one designed for TVD schemes. Hence the only difference in computation is that the TVD scheme reguires more elaborate dissipation terms or flux limiters for the explicit and implicit operators. In order to gain more insight into the TVD scheme for the complex flow problems, one objective of this

PAGE 21

10 dissertation is to examine how these various flux limiters affect the accuracy as well as efficiency of the numerical solutions. An extended derivation of the TVD schemes in threedimensional generalized curvilinear coordinates will be presented. The formula will then be applied to actual flow computations. A time-dependent axisymmetric Navier-Stokes code developed at the U.S. Army Ballistic Research Laboratory [26] has been modified to utilize more sophisticated numerical dissipation models based on the TVD schemes. Furthermore, a spatially varying time-step procedure [4] is also incorporated into this new version of the code to speed-up the convergence rate for steady-state flow computations. The transonic turbulent flow over the axisymmetric secant-ogive-cylinder-boattailed (SOCBT) projectile with sting at zero angle of attack is used as the test problem to assess the performance of the TVD schemes in terms of numerical accuracy, robustness and convergence rate. The transonic turbulent flow over an axisymmetric projectile problems at three different Mach numbers has been computed on many different grid sizes with both the Beam-Warming scheme and three variants of the TVD scheme. A systematic assessment will be presented to appraise the relative strength and weakness among these schemes.

PAGE 22

11 Numerical results will also be compared with the experimental data obtained by Kayser et al [27]. The remaining dissertation is organized as follows. The basic properties of hyperbolic conservation laws and TVD schemes are first reviewed in Chapter II. Then the detailed formulation of upwind and symmetric TVD schemes for nonlinear scalar hyperbolic eguations and system of hyperbolic eguations with constant coefficients are described in Chapter III. In Chapter IV, A brief survey of the mechanism of flux limiters, which play the central role in all TVD schemes, will be presented. In Chapter V, The derivation is presented which extends the TVD schemes in two-dimensions to three-dimensions with generalized curvilinear coordinates. In Chapter VI, a thin-layer Navier-Stokes solver along with some descriptions of the related physical modeling and numerical procedure aspects are presented. Results of numerical computation and detailed discussions are given in Chapter VII, which is followed by a summary and concluding remarks in Chapter VIII.

PAGE 23

CHAPTER II PRELIMINARIES In this chapter we give a summary of important notions in the theory of hyperbolic conservation laws and conservative finite difference methods. We also briefly review the theory of total variation-diminishing (TVD) finite difference schemes . 2.1 Hyperbolic Conservation Laws In this section, we will consider numerical approximations of the initial value problem for nonlinear hyperbolic systems of conservation laws: (2.1) U t + [f(u)] x =0 -co < X < co, u(x,0) = u (x) , where u(x,t) is a column vector of m unknowns, the flux function f(u) is a vector of m components, and u (x) is the initial data. Equation (2.1) can also be written in a matrix form as: (2.2) u t + A(u)u x = 0, A(u) = af/au. 12

PAGE 24

13 Equation (2.1) is called hyperbolic if the Jacobian matrix A(u) has real and distinct eigenvalues K l < K 2 < < * m Here k^, i=l, — ,m , are eigenvalues of the matrix A(u) . It is well known that the solution of Eq. (2.1) may develop discontinuities even when the initial data are smooth [11]. Because of this, we seek a weak solution which satisfies Eq.(2.1) in an integral sense, i.e., (2.3) J [u t x f (u) ]dx dt + [u (x)^]dx = 0, where (x,t) is a continuously differentiable test function. An immediate consequence of the integral relation Eq.(2.3) is that every piecewise continuous weak solution must satisfy a jump condition, which is known as the Rankine-Hugoniot condition in gas dynamics, across the discontinuity, i.e., (2.4) s [uj.] = [fi] i= i,..., m Here the brackets [ ] denote the jump across the discontinuity, and s denotes the speed of the discontinuity .

PAGE 25

14 It is also well known that weak solutions of hyperbolic conservation laws are not uniquely determined by their initial values. A simple example of such a discontinuous solution of Eq.(2.1) for a scalar conservation law with f (u) = 1/2 u 2 is ui(x,t) = The discontinuity is across the line x = t/2, which propagates with speed s = 1/2; across the discontinuity [u] = 1, [f] =1/2, so the jump condition of Eq. (2.4) is satisfied. The function 1 for X < t/2 for X > t/2 u 2 ( x.t, . { ; for x < t/2 for x > t/2 again is discontinuous solution across the line x = t/2, and [u] = -1, [f] = -1/2, so that the relation of Eq. (2.4) is satisfied also. In order to select the physically realizable solution among the weak solutions, some additional condition has to be satisfied. The criterion used is that the 'correct' solution should be the one obtained in the limit as e -* in the viscous equation (2.5) ut + f(u) x = eu XX'

PAGE 26

15 Lax [28] has defined an entropy inequality, with the aid of an entropy function U(u) for Eg. (2.1), which is reguired to satisfy the following two conditions: (i) U satisfies (2.6) U u f u = F u . where F is some function called entropy flux, (ii) U is a convex function of u. i.e., (2.7) U uu > 0. Multiplying Eg. (2.5) by U u , using the conditions of Eg. (2.6) and Eg. (2.7), also taking the limit as e + shows that the correct solution is the one satisfying (2.8) U(u) t + F(u) x < 0. If u is a weak solution with discontinuity travelling with speed s, then (2.9a) s« [U] < [F] . or (2.9b) s. (U R U L ) < (F R F L ).

PAGE 27

16 When applied to the gas dynamics, Eq. (2.8) or (2.9) is equivalent to the fact that entropy increases across the shock waves. For this reason we shall call Eq.(2.8) or Eq. (2.9) the entropy condition or entropy inequality. In the scalar case, Oleinik [29] gave the entropy condition across the shock as f(u) f(u L ) f(u) f(u R ) (2.10) < s < u u L u u R for all u between ul and ur. Notice that the role of the entropy condition is to attempt to distinguish those discontinuous solutions which are physically realizable from those which are not. Lax [28] then defined the kth field to be genuinely nonlinear if (2.11) (grad /cjj.ri * 0. where ri is the kth right eigenvector of A(u) . He also defined the field to be linearly degenerate if (2.12) (grad «i).ri = 0. For a genuinely nonlinear field, we call the discontinuity a shock or shock wave if for some index i

PAGE 28

17 (2.13) «i(u R ) < s < /cj.(u L ) . For linearly degenerate fields, we call the discontinuity a contact discontinuity if (2.14) *i(u R ) = s = /ci(u L ). In one spatial dimension, the inviscid equations of gas dynamics can be written in the conservative form as dq 3f(q) + = 0, at ax where q = (p, m, e ) T . is the vector of conservative variables, f = (m, (m 2 / P )+p, (e+pjm/p ) T is the flux vector, and m = pu. In the above equation p is the pressure, P is the density, u is the velocity, and the total energy per unit volume is denoted by e = p/(7~l) + m The eigenvalues of the Jacobian matrix A(q) = af/3q are

PAGE 29

18 «i = u c, «2 = U, K3 = u + c where c is the local speed of sound given by c = (7P/p)^. The corresponding eigenvectors are ri = (1, u-c, H-uc) T r 2 = (1, u, hn 2 ) T rj = (1, u+c, H+uc) T where H = (e+p )/ P = c 2 /( 7 -l)+ hu 2 is the enthalpy. The wave families associated with k± and k 3 can be rarefaction waves or shock waves. The wave associated with k 2 is a contact discontinuity. The contact discontinuity is due to the original discontinuity in the data; it is presented in the linear approximations to the equations. Therefore, the corresponding characteristic field is called linearly degenerate. The shock waves are the discontinuities due to the nonlinearities in the equations; the corresponding characteristic field is called genuinely nonlinear. We now consider an explicit (2k+l) -point finite difference scheme in conservation form approximating Eq. (2.1) n+1 n (2.15) ui = ui A[ h i+1/2 hi_ 1/2 ]

PAGE 30

19 where A = At/Ax h i+l/2 = h ( u i-k+l/ u i-k+2/---# u i+k) h i-l/2 = h ( u i-k/ ui_ k+1 ,..., Ui+k-x) and hi+i/2, the numerical flux, must be consistent with the physical flux in the sense that h(u,u,...,u) = f (u) . In 1960, Lax and Wendroff [1] showed that if v(x,t) is a solution of such a finite difference scheme in conservation form, and if v(x,t) converges almost everywhere to some function u(x,t) as ax and At tend to zero, then u(x,t) is a weak solution of Eq. (2.1). Conservation forms maintained by the finite difference schemes are very important since this ensures that the solutions generated by such schemes have the correct shock speed. Another important theoretical contribution was made by Harten et al., [6], who defined a class of monotone finite difference schemes. A scheme n+1 (2.16) ui = H(ui_ k , Ui_ k+1 ,..., Ui +k ) n = Ui -A [h(ui_ k+1 ,. .., Ui +k ) h(ui_ k , . . ., ui+j^!)]

PAGE 31

20 is said to be monotone if H is a monotonically increasing function of each of its arguments. All convergent solutions of such finite-difference monotone schemes will satisfy entropy conditions and hence are physically realizable solutions. Harten et al., [6] also showed that all monotone schemes in conservation form can only be first-order accurate. 2.2 Background of TVD Schemes Consider the scalar hyperbolic conservation law (2.17) u t + [f(u)] x = u t + a(u)u x = 0. A weak solution to this problem has the following monotonicity properties: (1) No new local extrema in u(x) may be created. (2) The value of a local minimum is nondecreasing and the value of a local maximum is nonincreasing. It follows from these monotonicity properties that total variation (TV) in x of u(x,t) does not increase in t; i.e., (2.18) TV(u(t 2 ) ) < TV (u(t 2 ) ), for all t 2 > t x . n Let ui be the numerical approximation of the solution of Eg. (2.17) at x =iAx and t = nAt, with Ax the spatial mesh size and At the time step. Consider a one-parameter family of five-point difference schemes in conservation form

PAGE 32

21 n+1 n+1 n+1 (2.19) ui + \6(h i+1/2 hi_ 1/2 ) n n n = ui A(l-0) (h i+ i/2 hi_ 1/2 ) n n where e is a parameter, A = At/Ax, hi +1/2 = h(ui_i n\, n n ' u i+l/ u i+2 ) / and function hi+i/ 2 is called the numerical flux function. Let n n+1 R i+l/2 (l-0)hi +1/2 + * h i+l/2 be another numerical flux function. Then Eq.(2.19) can be rewritten as n+1 n (2.20) Ui = ui A(fi i+1/2 fii_ 1/2 ) This new numerical flux now is a function of eight variables, B n n n n+1 n+1 n+1 n+1 R i+l/2 = h(ui_!, ui, ui +1 , u i+2 , ui.i, ui , u i+1 , ui +2 ), and is consistent with the physical flux f in the following sense R(u, u, u, u, u, u, u, u ) = f(u) To simplify the notation, Eq.(2.19) can be rewritten as (2.21a) L.u n+1 = R.u n

PAGE 33

22 or (2.21b) u n+1 = (L _1 .R).u n = S.u n where L and R are the finite difference operators: (2.22a) (L.u)i = ui + A0(h i+1/2 hi_ 1/2 ) (2.22b) (R-u)i = ui A(l-0) (hi +1 /2 -hi_ 1/2 ) and S = L _1 .R The total variation of a mesh function u n is defined to be (2.23) TV( U n) = I | ui +1 Ui| = | |A i+1/2 un| i=— oo i=— oo Here and throughout this studies, the general convention A i+l/2 (*) = (*)i+i ~ (*)i is used. The numerical scheme of Eq. (2.19) for an initialvalue problem Eq. (2.17) is said to be TVD if (2.24) TV(u n+1 ) = TV(S«u n ) < TV(u n ) In addition, a scheme is called monotonicity preserving if the finite difference operator S is monotonicity preserving; that is, if u n is a monotone mesh

PAGE 34

23 function, so is the S»u n . Following Harten [11], the TVD schemes are monotonicity preserving. This is the reason why a TVD scheme will not produce spurious oscillations. The following sufficient conditions for Eq. (2.19) to be a TVD scheme are due to Harten [11], (2.25a) TV(R.u n ) < TV(u n ) and (2.25b) TV(L.u n+1 ) > TV(u n+1 ) Assume that the numerical flux h in Eq. (2.19) is Lipschitz continuous and that Eq. (2.19) can be written as n+l _+ (2.26) Ui -A0(C i+1/2 A i+1/2 u Ci1/2 Ai1/2 U) n+1 n ._ + Ui + A (1-0) (Ci +1/2 A i+1 /2U C i _ 1 / 2 Ai_ 1/ 2U) n where C i+1/2 = c (u^, Ui , u i+1 , u i+2 ) and Ci_ 1/2 = C + ( u i-2/ u i-if ui, Ui +1 ) are some bounded functions. Then Harten [11] further showed that sufficient conditions for Eq. (2.24) are (a) if for all i ± ± (2.27a) C i+1/2 = A(l-0)C i+1/2 >

PAGE 35

24 + .+ (2.27b) C i+1/2 + C i+1/2 = A (1-0) (C i+1/2 + C i+1/2 ) < 1 and (b) if for all i _ + (2.28) -oo < c < X9 Ci +1 / 2 < for some finite C. This set of sufficient conditions are very useful in checking or constructing the TVD schemes which do not exhibit the spurious oscillations associated with the conventional second-order shock-capturing schemes.

PAGE 36

CHAPTER III UPWIND AND SYMMETRIC TVD SCHEMES In this chapter, the description of nonlinear, upwind and symmetric, explicit and implicit, second-order accurate, entropy satisfying schemes for hyperbolic conservation laws will be presented. The notion of TVD schemes was first introduced by Harten [5]. He developed a class of first-order TVD schemes that prevent the growth of the total variation of the solution, and therefore the possibility of spurious oscillations can be eliminated. He also devised second-order accurate TVD schemes by adding limited anti-diffusive flux to a first-order accurate entropy satisfying upwind scheme. His method is sometimes referred to as modified flux approach. A more complete description of the TVD scheme based on Harten "s modified flux approach will be presented in the following sections. 3.1 First-Orde r Accurate TVD Schemp Consider the one-dimensional scalar hyperbolic conservation law 25

PAGE 37

au af(u) (3.1) + 26 au + a(u) du • — at ax at ax and with the initial data (3.2) u(x,0) = u (x) where a(u) = df/du is the characteristic speed. As described in Chapter II, the discrete approximation of Eq. (3.1) by a three-point explicit conservative finite difference scheme can be written as n+1 n n n < 3 3 > u i = »1 " Mhi+i/2 " hi_ 1/2 ) where A = At/Ax. The numerical flux in Eq(3.3) can be defined as (3.4) h i+1/2 = ^[fi + f i+1 ^(a i+1/2 )A i+1/2 u] where f L = f( Ui ) and A i+1/2 u = u i+1 Ui/ and (3.5) ai +1 / 2 = i ( f i+l _ fi)/(Ai + i/ 2 u) if Ai +1 / 2 u * df / du i if A i+1/2 u = o, U = Ui In equation (3.4), rj, can be viewed as the coefficient of numerical viscosity, and it is a function of A and a i+1/2 .

PAGE 38

27 Since the solution of the difference scheme of Eq.(3.3) does not necessarily possess the TVD property, spurious oscillations (wiggles) may be produced, especially in the neighborhood of discontinuities (shock waves, contact discontinuities) . In order to eliminate such spurious oscillations of the numerical solution, Harten [5] has shown that if a finite difference scheme can be written in the form of n+l n + (3.6a) Ui = Ui AC i _ 1/2 Ai_ 1/2 U + AC i+1/2 A i+1/ 2U where + (3.6b) ACi_ 1/2 = A[ ai_ 1/2 + tf(ai_ 1/2 )]/2 (3.6c) AC i+1/2 = a[ -a i+1/2 + iKa i+1/2 )]/2 then the sufficient conditions for the numerical scheme to be TVD are (3.7) AC i+l/2 + > (3.8) AC i+l/2 > + ( 3 ' 9 > A (Ci+l/2 Ci+i/2) * lNote that Eq. (3.9) imposes a Courant-Friedrichs-Lewy (CFL)like condition on the difference scheme Eg. (3.3). Unlike the monotone schemes, TVD schemes do not automatically satisfy the entropy condition, and the schemes

PAGE 39

28 may converge to a non-physical (entropy condition violated) solution. Consequently, some mechanism may have to be added to the TVD schemes in order to enforce the selection of the physical solution. As shown in [11], a slight modification of the coefficient of the numerical viscosity term r m (3.10) *(z) = i 1*1 > (Z 2 + e 2 )/2 Z < £ can avoid violating the entropy condition. In Eq. (3.10) r/>(z) is an entropy correction to |z|, and i is a small positive parameter (see reference [11] for detailed formulation of e). It is well known that an explicit scheme subject to a CFL-condition is limited to use of a small time step. Therefore, in order to obtain steady-state solutions, it usually requires long computing time. This penalty is particulary severe for high Reynolds number Navier-Stokes calculations where the fine meshes needed to resolve small length scales presented in the flow field dictates the use of very small time steps. The implicit schemes, on the other hand, can usually adopt much larger time steps. In [18], Harten devised a one-parameter family of conservative schemes of the form (3.11) n+1 n+1 n+l Ui + A0(h i+1/2 hi_ 1/2 ) =

PAGE 40

29 n n n ui -A (1-0) (h i+1/2 hi_ 1/2 ) When 0=0, Eq. (3.11) is reduced to Eq. (3.3), an explicit method. When 6*0, Eq. (3.11) is an implicit scheme. For example, if 6 = 1/2, the time differencinq is the trapezoidal formula, and if 9 =1, the time differencinq is the backward Euler method. In qeneral Eq.(3.11) is first-order accurate in space. When 6 = 1/2, the scheme is second-order accurate in time. Equation (3.11) can be written as (3.12) L«u n+1 = R»u n where L and R are the finite-difference operators, defined as (3.13a) (L.u)i = ui + A0(h i+1/2 hi_ 1/2 ) (3.13b) (R-u)i = ui A(l-0)(h i+1/2 hi_ 1/2 ). A sufficient condition for Eq.(3.11) to be a TVD scheme is that (3.14a) TV(R.u n ) < TV(u n ) (3.14b) TV(L«u n ) > TV(u n+1 ).

PAGE 41

30 A sufficient condition for Eq. (3.14) is the CFL-like condition 1 (3.15) Uai+l/2l * AV>(ai +1/2 ) < i e It is noted that, according to Eq.(3.15), the backward Euler implicit scheme (0 = 1 in Eq. (3.11)) is unconditionally TVD, while the trapezoidal formula (8 = 1/2) is TVD under the CFLlike restriction of 2. The forward Euler explicit scheme {6 0) is TVD under the CFL restriction of 1. 3.2 Second-Order Accurate TVD Scheme In order to achieve a second-order accurate scheme which satisfies the TVD conditions, Harten modified the numerical flux fi of a first-order scheme by adding a limited antidiffusive flux gi to the numerical flux (3.16) t± t i .+ gi where g^ approximates au 1/2 Axy>(a) ax and thus cancels the first-order error. Therefore, Harten 's method is sometimes referred to as the modified-flux approach. The modified flux is chosen so that the scheme is second-order

PAGE 42

31 in smooth regions and first-order at points of extrema. After introducing the modified numerical flux, the fivepoint, second-order accurate backward Euler implicit TVD scheme (i.e., f = 1) of Eq. (3.11) can be written as n+1 _n+l _n+l n (3.17) ui + A(h i+1/2 hi_ 1/2 ) m ui and the modified numerical flux hi +1 / 2 can be defined as (3.18a) hi +1/2 = [fi + f i+1 \Hai +1/ 2 + 7i+l/2) M+l/2 u ]/2 (gi+l-9i)/(Ai + i/2U) if Ai +1 / 2 u * (3.18b) 7i+l/2 if a 1+1/2" = where £i = fi + gi It is necessary to limit gi to prevent the possibility of Ti+l/2 becoming unbounded. A simple choice of g^ is gi = S.max[0, min(ai +1/2 |Ai +1/2 u| , S.ai_ 1/2 Ai_ 1/2 u) ] where S = sign(Ai +1/2 u) and

PAGE 43

32 a (Z) = **V>(Z) > for steady-state calculation. One can rewrite g^ as gi = minmod(ai +1 /2Ai +1 / 2 u, <7i-i/2M-l/2 u ) where the minmod function is defined as follows: minmod (x,y) = sgn(x) •max{0, min[|x|, y«sgn(x)]} or minmod (x,y) = • if x and y have opposite signs x if |x| < |y| and x and y have the same sign y if |y| < |x| and x and y have the same sign The minmod function is used here mainly to avoid overshoots and undershoots in the vicinity of discontinuities. Therefore gi is the so-called 'flux limiter' function for the control of unwanted oscillations in the numerical scheme. After careful examination of Harten's modified-flux approach, Yee found that the numerical flux hi+i/ 2 used in Eq. (3.18a) is quite diffusive [12]. A slightly modified less diffusive numerical flux can be defined as follows: (3.19) h i+1/2 = [fi + f i+1 + * i+ i /2 ]

PAGE 44

33 where ^i+1/2
PAGE 45

34 at most first-order accurate. Neverthless, It is possible to extend the scalar scheme to system cases such that the resulting scheme is TVD for the " locally frozen" constant coefficient system. The idea behind this extension relies heavily on the Roe's approximate Riemann solver [30]. Consider a one dimensional hyperbolic system au 3E(U) (3.20) + = at ax or au au (3.21) + A(U) = at ax Here U and E(U) are column vectors of m components. Roe first linearized Eq.(3.21) to au . au + a = o at ax where A(U R ,U L ) is a constant matrix for any given Ul, U r which represent the local conditions. The initial data for this Riemann problem is U(x,0) = ' U L X < . U R X > The matrix A(U L ,U R ) satisfies the following properties which

PAGE 46

35 are the so-called Property U defined by Roe [15]. (1) A(U L ,U R ) is a linear mapping from vector space U to the vector space E. (2) As U L U R + U, A(U L ,U R ) * A(U) , where A(U) = 3E/3U. (3) For any U L , U R , A(U L ,U R ) (U L -U R ) = E L -E R . (4) The eigenvectors of A are linearly independent. The eigenvalues of A may be viewed as the wave speeds of the Riemann problem and the projection of (U L -U R ) onto the eigenvectors may be regarded as the jump between the intermediate states. Conditions (3) ensure that the exact solution is obtained if (U L ,U R ) satisfy the jump condition S(U L -U R ) = (E L -E R ) where S is the shock speed. Condition (3) shows that S is an eigenvalue of A. Since A is the mean value Jacobian , in one dimensional Euler eguations of gas dynamics, such a matrix is constructed by evaluating its Jacobian matrix with specially averaged cell interface values (denoted by subscript i+1/2) of density, velocity and total enthalpy: Pi+1/2 = (Pi) h (Pi+i) is ui+l^i+l)^ + ui(pi) h u i+l/2 = (Pi+l) h + (Pi)*

PAGE 47

36 H i+l/2 = (Pi+l) h + (Pi)' 5 from which the speed of sound can be calculated as 2 Ci+l/2 = <(7-l)[ H i+1/2 ^(u i+1/2 ) 2 ]} Roe's approximate Riemann solver is further generalized by Yee [25]. Her approach is sometimes referred to as the " local characteristic approach " . The procedure is to define at each point a local system of characteristic variables W and to obtain a system of uncoupled scalar equations. After that, one can apply TVD algorithm of a scalar hyperbolic conservation law to each of the locally defined (frozen coefficients) characteristic variables of Eq. (3.20) as follows: (3.22) + A = 0, aw + A aw at ax w = R' -1 U A = R _1 UR = diag(a 1 ) where diag(a 1 ) denotes a diagonal matrix with diagonal elements a 1 which are the eigenvalues of the matrix A(U) , R denotes the matrix whose columns are eigenvectors of A, and

PAGE 48

37 R -1 is the inverse of R. A two-dimensional hyperbolic system of conservation laws is considered here to demonstrate the construction of the scheme , (3.23) U t + E(U) X + F(U) y = Here U, E(U) and F(U) are column vectors of m components. Let A = aE/aU and B = aF/au be the Jacobian matrices. Let the eigenvalues of A be (a x , a x , ,a x ) and the eigenvalues of B be (a y , a y , ,a y ) . Denote R x and R y as the matrices whose columns are the eigenvectors of A and B, and -1 -1 denote R x and R y as the inverses of R x and R y . Let the grid spacing be denoted as Ax and Ay such that xi = iAx and yj = I)Ay. Let ai +1 / 2 , Ri+i/2> R i+l/2 denote guantities a x , R x , R x evaluted at Ui +1 / 2 ,j respectively. Similarly, let aj +1 / 2 / —1 p _i R j+l/2» R j+l/2 denote the guantities a y , R y , R y evaluated at Ui^j+i/2 respectively. Define (3.24a) ai+l/2 = R i+l/2 ( u i+l,j " u i,j) as the difference of the characteristic variables in the local x-direction, and define (3.24b) "j+l/2 = R j+l/2 ( u i,j+l " Ui,j)

PAGE 49

38 as the difference of the characteristic variables in the local y-direction. With the above notation, a one-parameter family of TVD schemes (3. 12) in two dimensions can be written as n+1 x _n+l _n+l (3.25a) U^j + A 0(E i+1/2/ j Ei_ 1/2 ,j) y _n+l _n+l + A *(Fi f j+i/ 2 " Fi^j.i/2) n x _n .n = "i,j A (1 -0(Ei+i/ 2/ j Ei_ 1/2 ,j) y .n _n A (1*)(Fi f j+i/ 2 Fi^-i/2) with A x = At/Ax, Ay = At/Ay. Then the numerical flux function ^x+1/2,-} for hoth the upwind and symmetric TVD schemes can be expressed as (3.25b) E i+1/2/ j = %£«i f j + E i+lf j + Ri+i/ 2 *i+l/2]Similarly, one can define the numerical flux F^ j+i/23.4 Upwind TVD Schemes The elements of $i+i/2/ needed in Eq.(3.19) for a secondorder upwind TVD scheme, originally developed by Harten and later modified and generalized by Yee [12], are denoted by (^i+1/2) 1 £ " 1 , ,m and defined as follows:

PAGE 50

(3.26a) ( *i + i /2 ) U = 39 I I ( a i+l/2) (9i +
PAGE 51

40 Equation (3. 27b) and Eq. (3.27c) are suqgested by Roe[31] and van Leer [32] respectively. 3.5 Symmetric TVD Schemes In 1984, Davis [23] showed that Roe-Sweby's second-order TVD scheme can be interpreted as a Lax-Wendroff scheme plus an upwind weighted numerical dissipation term. A specific flux limiter then can be constructed by satisfying the sufficient conditions of TVD properties mentioned previously. Also the requirement of upwind weighting can be removed. Thus TVD schemes are not necessarily upwind any more, these non-upwind TVD schemes are referred to as symmetric TVD schemes. Shortly after that, Roe [24] reformulated Davis's approach in such a way that the formulation is easier to analyze. Flexible choices of flux limiters are also provided to enhance the shock resolution. Subsequently Yee [25] further generalized the works of both Davis and Roe to a class of symmetric TVD schemes. The formulation of Roe-Davis can be viewed as one of the variants of the symmetric TVD schemes proposed by Yee. The advantages of Yee's approach are that (1) stiff problems can be treated properly by using implicit schemes, (2) numerical boundary conditions are not so sensitive in comparison with upwind TVD schemes, (3) computational complexity can be further reduced.

PAGE 52

41 I S The elements of the ($1+1/2) denoted by (^i+1/2) for a general second-order symmetric TVD scheme [25] are I i i J! (3.28a) (^i+1/2) -V>(ai+l/2) ["i+1/2 " Qi+1/2 I« Examples of the 'limiter 1 function Qi+1/2 can be expressed as i 1 J I J (3.28b) Qi+1/2 = minmod(ai_i/ 2 ,ai + i/2) + minmod(ai + i/ 2 ,ai + 3/ 2 ) I " a i+l/2 I 2 I I (3.28c) Qi+1/2 = minmod(ai_i/ 2 , a i+l/2f a i+3/2) 1 a a a (3.28d) Qi+1/2 = minmod [ 2ai_i /2 , 2ai+ 1/2 , 2ai +3 / 2 , I s. l/2(a-i-i/2 + ai+3/2)] i i where ai+1/2 and ^(ai + i/ 2 ) are defined as before. In general, the limiter of Eq. (3.28d) is more compressive than the other two. Since Eq(3.25a) is a highly nonlinear implicit scheme, in order to gain computational efficiency, local linearizations are applied to the nonlinear terms (denoted by n+1 time level) of Eq. (3.25a). There are two linearized forms of Eq. (3.25a) proposed by Yee. The first method will maintain the

PAGE 53

42 conservative property, but the resulting scheme may no longer be unconditionally TVD. This method is referred to as the linearized conservative implicit (LCI) form. The second method will destroy the conservative property but preserve the unconditionally TVD property. The latter method is referred to as the linearized non-conservative implicit (LNI) form. However, according to Yee and Harten [12], It appears that the LCI form had better convergence rate than the LNI form. Thus the only LCI form is adopted here. 3.6 Linearized Conservative Implicit (LCI) Form The nonlinear terms of Eg. (3.25a) are linearized in time about U n by a Taylor series such that E n+1 = E n + A n (u n+1 _ ^ + (At 2 ). Next we can locally linearize the coefficient of Ai +1 / 2 U by dropping the time level from (n+1) to n, then obtain the LCI form in delta formulation as XX XX (3.29a) [I + A *H i+1/2 ,j A *Hi_ 1/2 ,j y y + Ay^H i/j + 1/2 Ay
PAGE 54

43 _n _n ^ y [Fi,j+l/2 " ^1^-1/2 ] where x x (3.29b) H i+1/2/ j = ^[A i+1/j + Oi+l/ 2 ,j ] n Y Y (3.29C) H i/j+1/2 = ^tBi,j+i + ni,j+l/2 ] n Here A and B are the Jacobians of the fluxes E and F, and (3.30a) Ai+i/ 2 ,j = ( R x diag[ p*iKa* + y*) ]Rx*) 1+1/2*1+1/2 y I i l -1 (3.30b) ni^+i/2 = ( Ry diag[ *(a + 7 ) ]R y )j+i/2Aj+i/ 2 where (3.30c) 4+1/2 = (si + gi+l)/« 1+1/2 » i i, 2 a £j+l/2 = (9j + gj+l)/"j+l/2 For steady-state application, one way to simplify Eq.(3.30) is to use a spatially first-order implicit operator; e.g., by redefining Eq(3.30a) and Eq. (3.30b) as

PAGE 55

44 (3.31a) «i+i/2,j = (R x diag[Vfa 1 )] Rx 1 ) i+i/ 2 Ai + i/ 2 (3.31b) n?,j+i/2 = (Ry diag[^(a 1 )] Ry 1 ) j+i/ 2 Aj + i/ 2 x The computation can be reduced even more if fii+i/2,j and v ... fli+j+1/2 are simplified to diagonal matrices. For example, one can redefine Eg (3. 31) as (3.32a) 0i+i/2,j = ( diag[ max ^(a 1 ) ] ) i+i/2Ai+i/ 2 y . (3.32b) "i,j+i/2 = ( diag[ max ^(a 1 ) ] ) 1+1/2* j+1/2* I Eguation(3.29) together with Eg (3. 32) is referred to as the linearized conservative diagonal form. 3.7 Alternating Direction Implicit (API) Form The integration of the two-dimensional implicit operators, such as Eg. (3.29a), is still very expensive. One can simplify the solution process by introducing an approximate factorization of the two-dimensional operator into two one-dimensional operators. This ADI form of Eg. (3.29a) reads

PAGE 56

45 X X (3.33a) [ I +A x 0H i+1/2 ,j A x 0Hi_ 1/2 ,j ] D* = A x [E i+1/2 ,j Ei-i/ 2 ,j] " AY tFi,j+i/ 2 " Fi,j-l/2] y y (3.33b) [ I + Ay^Hi^+i/2 ^0^1,^-1/21 D = D* (3.33c) U n+1 = U n + D The solution algorithm now consists of two onedimensional sweeps, one in xand one in y-direction. Each sweep requires the solution of a linear system involving a block tridiagonal matrix which can be solved very efficiently by block LUD (lower-upper decomposition) solver.

PAGE 57

CHAPTER IV FLUX LIMITERS As already mentioned, first-order accurate schemes, in general, suffer from numerical diffusion, white the classical second-order accurate schemes such as the Lax-Wendroff scheme [1], exhibit the spurious oscillations near discontinuities of the solution. In order to overcome this problem, One can derive the flux limiters where the anti-diffusive fluxes are limited so as to prevent these oscillations from occurring. In this chapter, the mechanisms of a class of flux limiters are analyzed. Some specific flux limiters which have been proposed by various authors are also examined. Consider the initial value problem for a scalar hyperbolic conservation law 3u 3f(u) (4.1) + = 0, at ax and initial value u(x,0) = u (x) . As before, the general three-point explicit finite difference in conservation form which approximates Eq.(4.1) can be 46

PAGE 58

47 written as n+1 n n n (4.2) ui = ui A(h i+1/2 hi_ 1/2 ) where satisfying h i+l/2 = h(ui +1 ,ui) , (4.3) h(u,u) = f(u), that is, the numerical flux h(u,u) should consist with flux function f (u) . If the scheme of Eq. (4.2) can be written in the form n+1 n n n (4.4) Ui = Ui C i _ 1/2 Ai_ 1/ 2U + Di +1 /2Ai +1 /2U , where Ci_i/ 2 and Di+i/ 2 are functions of u n , then it can be shown [5] that sufficient conditions for the scheme to be TVD are the satisfaction of following inequalities: (4.5) C i+1 / 2 > 0, (4.6) D i+1 / 2 > 0, (4.7) C i+1/2 + D i+1/2 < 1. For easy presentation, the following scalar linear equation is considered

PAGE 59

48 3u an (4.8) + a =o a = const. > 0. at ax In [34], Sweby has shown how the second-order Lax-Wendroff scheme could be constructed by adding an anti-diffusive flux to the first-order upwind scheme, that is n+1 n n n (4.9 ) ui = Ui ^Ai_ 1/2 u a_[ H v(l-v)A i+1/2 u ], where CFL number w Aa, and a_ is the backward difference operator, i.e., a_u = ui *£«.].. It is well known that the Lax-Wendroff scheme produces spurious oscillations near the discontinuities of the solution, while the first-order upwind scheme does not. By limiting the amount of the anti-diffusive flux in Eq. (4.9), it is possible to obtain the second-order accurate TVD scheme. The concept of this anti-diffusive flux is much the same as the flux corrected transport (FCT) approach proposed by Boris and Book [9], and Zalesak [10]. The main differences between the present flux limiters and FCT are that the former are single step schemes, the limiter being functions of data u n , whereas FCT is a two-step procedure. The anti-diffusive flux (the second term on the right hand side of Eg. (4.9)) is modified by introducing a limiting function ^,

PAGE 60

49 (4.10) -A_[ h i v (l-v)Ai +1 / 2 U ] The ^i is chosen to be a function of the ratio of consecutive cell gradients of A i-1/2 U (4.11) a = tf(ri) where ri = A i+1/2 U This parameter r^ plays a central role in most TVD schemes, The resulting scheme of Eq.(4.9) can be rewritten as n+1 n n (4.12) Ui = Ui u{l + %(l-v)[*(ri)/ri 4>{ri-i)-\)Li1/2 u • This is a scheme in the form of Eg. (4.4) with (4.13) Ci_ 1/2 = * (r) > for r > and that (f> (r) =0 for r < 0. Under these additional restrictions, the

PAGE 61

50 bound of Eq. (4.15) becomes (4.16) < {r)/r < 2; < (r) < 2. Note that if (r) = 1, Eq.(4.12) reduces to the central difference Lax-Wendroff scheme, and if 4> (r) = r, Eq. (4.12) reduces to the second-order upwind scheme proposed by Warming and Beam [35] . The region defined by Eq. (4.16) is shown in Figure 4.1. Another constraint for the scheme to be second-order accurate is imposed by requiring that the amount of limited flux added to the first upwind scheme be a weighted averaged of the second-order centered Lax-Wendroff scheme and the second-order upwind scheme of Warming and Beam. The corresponding region with this extra constraint is shown in Figure 4.2. In [32], van Leer averaged two second-order schemes to produce a second-order TVD scheme. He combined a monotonicity preserving but non-conservative version of the Lax-Wendroff scheme and the second-order upwing scheme of Warming and Beam to give a conservative monotonicity preserving version of Fromm's scheme [33]. The parameter he used as a " smoothness monitor " is defined by

PAGE 62

51 Warming-Beam r)=r Figure 4.1 TVD region

PAGE 63

52 Figure 4.2 Second order TVD region

PAGE 64

53 *i+l/2 A i+1/2 U Ai_ i-1/2 u which is the reciprocal of the riAfter reformulating the van Leer's scheme, the equivalent flux limiter is (4.17a) *(r) = r + r |r| + 1 2r (4.17b) +(r} = r > 1 + r r < which is a smooth curve and lies in the second-order TVD region of Figure 4.3. In [17], Roe proposed a second-order accurate monotonicity preserving scheme using an incremental (or fluctuating) approach. This involved the allocation of cell flux difference to update the values of u at the ends of the cells to obtain a first-order upwind scheme at the (n+l)At time level. Then a fraction of the increment is transferred across the cell against the direction of flow to achieve second order.

PAGE 65

54 i 2 Figure 4.3 Van Leer's limiter Eq. (4.17)

PAGE 66

55 The resulting scheme can be proved to be monotonic under certain restrictions on (r) . The most commonly used limiter is the " minmod" limiter which corresponds to the lower boundary of the TVD region in Figure 4.4, and which may be written as (4.18a) 4>{r) = max [ 0, min (l,r) ], or (4.18b) *(r) = r < [ min(l, r) r > Furthermore, Roe [31] has experimented with various transfer functions and proposed his more compressive "superbee" transfer function, which when written as a limiter is expressed as (4.19a) or (r) = max [0, min(2r,l), min(2,r)], (4.19b) *(r) P"« min(2r,l) min(2,r) r < < r < 1 1 < r which corresponds to the upper boundary of TVD region (see

PAGE 67

56 Figure 4.4 Minmod liniter, Eq. (4.18)

PAGE 68

57 Figure 4.5) In general, Roe's "superbee" limiter is the most compressive among the above three <$> (r) functions. Limiter (4.17b), due to van Leer, is found to produce a slightly better shock resolution than minmod limiter function Eg. (4.18a). m certain cases [31], "superbee" limiter is too compressive and might produce a non-physical solution. Therefore using different limiters for different characteristic fields could be a good compromise. For the one-dimensional Euler eguation of gas dynamics, the characteristic fields consist of two nonlinear fields u±c and one linear field u. One can further improve the sharpness of the contact discontinuities by using a slight diffusive limiter such as Eq. (4.17) for nonlinear fields and using a more compressive limiter such as Eq. (4.19) for linear fields. The above idea can be extended to a nonlinear scalar equation as well as to a multi-dimensional system of conservation laws. In 1984, Davis [23] interpreted the upwind TVD finite difference schemes analyzed by Sweby [34] as Lax-Wendroff schemes plus an upwind artificial dissipation term. Since the determination of upwind directions for hyperbolic systems is a complicated procedure, he modified this artificial dissipation term in such a way to make this term independent of the wave directions. The result is a dissipation term

PAGE 69

58 4>(r) 32Roe's superbee limiter 1 2 3 Figure 4.5 Roe's superbee limiter, Eq(4.19) -

PAGE 70

59 which is based on TVD constraints and which does not contain any upwind weighting. Shortly thereafter, Davis •s work was reformulated by Roe in a way which is easier to analyze. The analysis exhibited a class of TVD schemes which was not observed by Davis. Roe's scheme can be written as n+1 n (4.20) Ui = Ui hv (l+j/)Ai_ 1 / 2 U -hv(l-v)Ai +1 / 2 u %|*j (1-1* I) (i-Qi-i/2)Ai-i/ 2 u h\u\ (1-M) (l-Qi+l/2)Ai+l/2U. Here the first line represents the original Lax-Wendrof f scheme, and the other terms represent conservative dissipation terms. The function Qi+1/2 depended on three consecutive gradients Ai_;iy 2 u, Ai +1 / 2 u and Ai +3 / 2 u is of the form + (4.21a) Qi+i/ 2 = Q( r i+1/2 , r i+1/2 ), where 7 A i-V2 u + A i+3/2 u (4.21b) r i+1/2 = ; r i+1/2 A i+1/2 U A i+1/2 U Roe further showed that Q can be chosen in such a way to ensure that Eg. (4.20) is TVD, That is Eg. (4.20) can be of the form of Eg (4. 7) with

PAGE 71

60 (4.22) Ci_ 1/2 = i/[l-^(l-i/)Qi_i/ 2 +^(l-^)Qi+l/2/ri+l/2 1# D i+l/2 = 0. if a > 0, and the case for a < follows in a similar way, i.e. , (4.23) Ci_ 1/2 = 0, Di+l/2 1*1 [1 %.(1 IH) Qi+1/2 + H(l M)Qi-l/2/*i-l/2 ]• If one assumes that both Q and Q/r are always positive, then a set of sufficient conditions for Eq. (4.20) to be TVD are 2 (4.24a) Qi+l/2 < 1 -\u\ 2 (4.24b) Qi+1/2/ ri +1 /2 < M + 2 (4.24C) Qi+1/2/ ri +1/2 < M Based on these bounds, some flux limiters Q can be defined as (4.25a) Q(r , r + ) = minmod(l, r") + minmod(l,r + ) -1,

PAGE 72

61 (4.25b) Q(r~, r + ) = minmod(l, r~, r + ) , (4.25c) Q(r _ , r + ) = minmod [2, 2r", 2r + , h (r~ + r + ) ] . The graphical representations of these three limiters for symmetric TVD schemes are shown in Figure 4.6(a,b,c).

PAGE 73

62 o r"! 1 l r + -l r"+r + -l r + -1 1 r"-l (a) Q(r~,r ) = minmod(l,r~) + minmod(l,r ) 1 i (b) Q(r~,r + ) rainmod(l,r~,r ) (c) Q(r ,r ) = minmod(2,2r~,2r + ,l/2(r~+r + )) Figure 4.6 Graphical representation of. three limiters, Eq(4.6 a,b,c)

PAGE 74

CHAPTER V ON THREE-DIMENSIONAL SYMMETRIC AND UPWIND TVD SCHEMES IN GENERALIZED COORDINATES An original derivation of the implicit, second-order accurate symmetric and upwind TVD schemes for threedimensional hyperbolic systems of conservation laws in generalized curvilinear coordinates will be presented in detail in this chapter. 5.1 The 3-D Euler Equation in Generalized Coordinates The inviscid three-dimensional conservation equations of mass, momentum, and energy can be represented in a fluxvector form that is convenient for numerical simulations as (5.1a) aq + 3E(q) + 3F(q) ay + 3G(q) dz at ax o, where 63

PAGE 75

64 (5.1b) q = p pU. pW pW e E = p\l pu 2 + p puv , F = pUW u(e+p) pv pVU pv 2 + p pVW v(e+p) G = p\l pWU pWV p\t 2 + p w(e+p) The primitive variables of Eq. (5.1b) are the density p, the velocity components u, v and w, and the pressure p. The total energy per unit volume e is related to p by the equation of state for a perfect gas (5.1c) p = (7-1) [e hp (u 2 + v2 + w 2 )] where 7 is the ratio of specific heats. A generalized coordinate transformation of the form (5.2) r = t, f = £ (t, x, y, z) is the logitudinal coordinate v = ri(t, x, y, z) is the circumferential coordinate f '» t (t, x, y, z) is the near normal coordinate which maintains the strong conservation law form of Eq. (5.1a) is expressed as (5.3) aq dr + at(g) 3£ *F(q) dr, + aft(q) ar = o,

PAGE 76

where A q I I g 65 q/J, (It q + C x E + ^ y F + | z G)/J, (»?t q + »/x E + »?y F + f| G)/J, (ft q + fx e + r y f + f z g)/j, and 3(x, y, z) J = , the Jacobian of transformation. dtt, r# D The strong conservative form of the transformed equations is maintained in order to admit shock capturing. The contravariant velocity components along £-, r/-, and f -coordinates are defined as follows (5.4a) U=?t + ^U+^yV+^W, V = r/ t + »? X u + »?V V + ^Z w / w = ft + Tx u + fv v + fz w Here E, F and G of Eg. (5.3) becomes (5.5) E = pU pV pW 1 puU+ £ x p 1 pUV+»? x p 1 puW+f x p J pVU+ £yP pWU+ £ z p (e+p)U-| t p , F = J pW+r/yP pWV+r/zP (e+p)V-»7 t p ,G = J pVW+fyP pwW+r z p (e+p)W-f t p

PAGE 77

66 Let A = 3E/aq, B = 8F/dq and C = dG/dq; then the Jacobians A, B and C of ft, £ and 6 can be written as (5.6a) A = £ t i + £ x a + £ y b + £ z c, (5.6b) B = vt I + ^ x A + r) Y B + nz C, (5.6c) C = f t I + rx A + f y B + f z C. The detailed derivation of three-dimensional compressible inviscid gas dynamic equationss in generalized curvilinear coordinates can be found in Pulliam and Steger [36]. Let c = (tP/p)^ be the local speed of sound, The eigenvalues of A are (5.7a) (a^, a^, a c , a^ , a ? ) = (u, U, U, U+k^c, U-k^c) where k^ = (£ x 2 + £ y 2 + ^2^. The eigenvalues of B are (5.7b) (ar,, ar,* a% , an' a?) = (V, V, V, V+k^c, V-k^c) where k„ = (r, x 2 + „ y 2 + r, z 2 )^. The eigenvalues of t are (5.7c) (a r , a r , a r , a*, af ) = (w, W, W, W+k f c, W-k f c) where k r = (rj£ 2 + fy 2 + tm *)h m

PAGE 78

67 12 3 4 5 Furthermore, let R^ = (R^, R^ , R^ , R^ , R^ ) , R„ = (rJ, R 2 ,, rJ, rJ, rJ), and R f = (rJ, R 2 , r£, R*, R^) be the matrices whose columns are eigenvectors of k, & and 6. —1 —1 —1 Let R^ 7 R^ , and R f be the inverses of R^ , R^ , and R f . The forms of the matrices, R K with k = £, r; , or f and their inverse can be written as (5.8) K *x k y k z 1 1 k x u ky\l-k z p k z U+kyp u+k x c u-k x c k*v+k z p /CyV k z V-k x p V+kyC V-kyC k x w-kyp kyW+k x p ic z w W+k z C w-k z c k x q 2 /2+ P (k Z V-icyW) k y q 2 /2+ p (k x w-k z u) k z q 2 /2+ P (ic y u-k x v) H+C0 E-c'e (5.9) -1 R* = K x( 1 ~ lD l)~ «X b 2 u K z /p + ~Ky/p+ -«x b 2 (« Z V-KyW)/p " X b 2 v « x b 2 w /Cy(l-bl)~K z /p + Kyb 2 V *x/p + -Kyb 2 («; x W-/c z U)/p /Cyb 2 U Kyb 2 W Kz( 1_b l)" Ky/p + -/c x /p + /c z b 2 w -«z b 2 (»CyU-/C X V)/p *Z b 2 u « z b 2 v (b 1 -s/c)/2 (k x /C(Ky/C(/c z /Cb 2 /2 b 2 u)/2 b 2 v)/2 b 2 w)/2 (b 1 +fl/c)/2 -(/C X /C+ -(/Cy/C+ -(/c z /c+ b 2 /2 b 2 u)/2 b 2 v)/2 b 2 w)/2

PAGE 79

68 where "x = K x/( K x 2 + «y 2 + «z 2 )^/ Ky = «y/(« x 2 + (Cy 2 + k z 2 )' 2 , «z = «z/( /c x 2 + *y 2 + «z 2 )^/ 5 = K X U + /CyV + /C Z W, q2 = u 2 + v 2 + w 2^ H = c 2 /(7-D + q 2 /2, bl = *>2 (q 2 /2), t>2 = (7-l)/c, C = (7P/p)*, Let the grid spacings be denoted by A£ , Arj and Af such that $ = iA£ , r/ = jA»? and f = kAf . Denote qi+i/2,j,k as some symmetric average of qi+i,-),* and qi,j,k ( for example, Roe's average). Let a{ +1 / 2 , Ri+i/2/ R i+l/2 denote the quantities of a^ , R^ , R ? related to A evaluated at qi+i/2,j,k I -1 Similarly, let aj +1 / 2 / R j+l/2f R j+l/2 denote the quantities of a^ , R,j , R^ related to I evaluated at qi,j+i/2,k and a k+l/2' R k+l/2/ R k+l/2 denote the quantities of a^ , R r , rJT 1 related to £ evaluated at qi,j,k+l/20ne can define

PAGE 80

69 (5.10a) ai+i/2 = R i+l/2 (<3i+l,j,k gi,j,k) ^( J i+l,j,k J i,j,k) as the component of (qi+i,j,k " <3i, j ,k) (omitting the j and k indices) in the locally i-th characteristic £ -direction. Denote -1 1 R j+l/2 (qi,j+l,k " <3i,j,k) (5.10b) a j+1/2 *-= • — , *(Ji,j+i,k Ji,j,k) as the component of (qi,j+i,k " 3i,j,k) (omitting the i and k indices) in the locally i-th characteristic r? -direction. And denote tn «**, J R k+i/2 (qi,j,k+i qj y j f k) (5.10c) «k+l/2 = • *(Ji,j,k+l Ji,j,k) as the component of (qi,j,k+l qi,j,k) (omitting the i and j indices) in the locally i-th characteristic f -direction. I The explicit expression for vector <*i+i/2 of Eq. (5.10a) is given by

PAGE 81

(5.11) 70 1 "i+1/2 <*i+l/2 <*i+l/2 4 a i+l/2 5 a i+l/2 -/CyAi +1 / 2 W+k z Ai +1 / 2 V+ 2 . 2 [ c i+l/2 A i+l/2^ _A i+l/2P]«x/ c i+l/2 *x A i+l/2 w -*z A i+l/2 u+ 2 . 2 [ c i+l/2 A i+l/2P~ A i+l/2P]«y/ c i+l/2 -"x A i+l/2 v+K y A i+l/2 u+ 2 . 2 [ c i+l/2 A i+l/2^A i+l/2P]«z/ c i+l/2 {Pi+l/2 c i+l/2( K x A i+l/2 u+K y A i+i/2 v+K z A i+l/2 w ) + 2 2 [ci+l/2 A i+l/2P -A i+l/2P] >/ c i+l/2 (/'i+l/2 c i+l/2(*x A i+l/2 u+ «y A i+l/2 v+ «z A i+l/2 w ) + 2 2 [ c i+l/2 A i+l/2^ _A i+l/2P] )/ c i+l/2 where (5.12a) A i+1/2P Pi+l,j,k " Pi,j,k, (5.12b) A i+1/2 U = ui +1 ,j, k ui^^k, (5.12c) >i+l/2 v = Vi +1/ j /k v i( j /kl

PAGE 82

71 (5.12d) Ai +1/2 W = W i+lf j /k Wi^j^k, (5.12e) M+1/2P = Pi+i,j,k " Pi,j,k* and Roe's form of the averaging in the £ -direction is: (5.13a) Pi+1/2 = (*>i+i,j,k)* (Pi,j,-k) H i is i^i „. X * Ui+1 '3' k + ?*>*}** (5.13b) Ui +1 / 2 ,j f k " with (5.13c) v i+1/2 ,j /k = (5.13d) W i+1/2/ j /k = X* + 1 X* v i+1/j , k + v i/j/k X* + 1 x * wj+l,j,k + w ifjyk X* + 1 (5.13e) Hi +1/2/ j /k = X* + 1 (5.13f) C i+1/2/ j, k = (7 1) [ H i+1/2f j /k h> (Ui+l,j /k + v? +lfj/k + wf +1/ j /k )] (5.13g) X* = (p i+ i,j /k /Pi,!,*)*, (5.13h) H = yp/p(y-l) + h(u 2 + v 2 + w 2 )

PAGE 83

72 It is noted that here H represents the total enthalpy and should not be confused with the numerical flux function H. Therefore, to use Roe's averaging for the £differencing, all one has to do is to compute n±+i/2 j k» v i+l/2,j,k/ w i+l/2,j,k and c i+l/2,j,k in Egs.(5.7), (5.8), (5.9), (5.10) and (5.11). Similarly Roe's averaging can be applied for the »?and f -directions. 5.2 Algorithm of 3-D TVD schemes in generalized Coordinates With the above notation, a one-parameter family of TVD schemes can be written as (5.14) qi/j,k + * * [Ei+i/ 2 ,j,k ~ Ei_ 1/2 j f k] + A^[F n ^ +1/2/k F n ^. 1/2 , k ] + ^n~Gi + ,],-k+l/2 G n ^,k-l/2] q?,j,k + (1 " 'H«!+l/2,j,lC " E n -i /2/ j /k 3 + (1 0[F n /j+V2/k ~F^ t j1/2l -k] + (1 9)[Gi t j fK+1/2 Gi f j /k _ 1/2 ]

PAGE 84

73 where S is a parameter, A* = Ar/A^, \ n = At/Atj and A^= Ar/Af. A particular form of the numerical flux function Ei+1/2 j k f° r both the upwind and symmetric TVD schemes can be expressed as (5.15) E i+1 / 2 ,j,k = %( &i t 1 t lt * ^i+l,j,k Ri+1/2 *i+l/2]Similarly, one can define the numerical fluxes F^ j+1/2 k and G i,j,k+l/2Upwind TVD Schemes I The elements of the $1+1/2 denoted by (^i+i/2) U for a second-order upwind TVD scheme, originally developed by Harten [5], and later modified and generalized by Yee [12], are I I J J (5.16 ) (^i+i/2) u = ^(ai+1/2) (gi+i gi) a i. 1 V>(ai + i/ 2 + 7i+i/2)<*i+l/2For steady-state calculation the function a(z) = l/2V»(z), and the functions ip and 7 are defined as before, i.e., (5.17) rf>(z) = and 1*1 |z| > . (Z 2 + £ 2 )/2e |z| < e

PAGE 85

(5.18)7i+l/2 = *(ai+l/2) < 74 1 a a (9i+l " gi)/«i+l/2 "i+1/2 *° ai +l/2 =0 Examples of the "limiter" functions g^ can be expressed as £ in. (5.19a) gi = minmod (a±1 / 2 i «i+l/2) / a a sl i. i (5.19b) gi = (ai +1/2 «i-l/2 + l<*i+l/2 »i-l/2l)/ i i («i+l/2 + a i-l/2) / (5.19c) gi = S. max[0, min(2 l«i +1/2 I , S'al-i/2), with min( |ai +1 / 2 | , 2S»ai_ 1 / 2 )], S = sign (ai+1/2) Again, the minmod function of Eq. (5.19a) is defined as follows minmod(x,y) = sgn(x)«max{ 0,min[ |x|, y«sgn(x)]} Symmetric TVD Schemes The elements of $1+1/2, denoted by (^i+i/2) S for a general second-order symmetric TVD scheme [25], are

PAGE 86

75 (5.20) (^i+i/2) S " " ^(ai+1/2) t «i+l/2 " Qi+1/2 ] A Examples of the "limiter" functions Qi+1/2 can be expressed as (5.21a) Qi+1/2 = minmod(ai_ 1 / 2 f a i+l/2) + minmod(ai +1 / 2 , <*i+i/2) "i+1/2 > A i a a sl (5.21b) Qi+1/2 = minmod(ai_ 1 / 2 , <*i+i/2/ "i+1/2) > k 4 J (5.21c) Qi+1/2 minmod[2ai_ 1/2 , 2ai +1/2/ I i i 2a i+l/2 / ^( a i-l/2 + "i+1/2)] 5.3 A Conservative Linearized API Form for Steady-State Application A conservative linearized ADI form of Eq. (5.14) for steady-state application can be written in the generalized coordinates as i i i i (5.22a) [I + A H i+1/2 ,j, k A 9 Hi_ 1/2 ,j /k ] D** i ~n _n " A [Ei +1 / 2/ j, k Ei_ 1/2/ j /k ]

PAGE 87

76 ~ A [ F i,j+l/2,k F if j-i/ 2/ k] f -n _n ~ A [ G i,j,k+l/2 ~ Gi / j /k _ 1 / 2]/ (5.22b) [I + A | H if j+i /2 ,k " A I IlJ, j-l/2 /k ] D* = D**, (5.22C) [I + A 9 H if j /k -fl/2 " A I H if j, k -i/ 2 ] D = D*, (5.22d) qn+1 = -n + D/ where (5.22e) H i+i/2,j,k = *f*l+i,j,k + ni+i/ 2 ,j,k] n , 9 (5.22f) H i/j+1/2/k = %t«J,j+l,Jc + "i,j+l/2,k] n , r r (5.22g) H i/j/k+1/2 = ^[6 i>j/k+1 + Oi,j, k+1 / 2 ] n , and € i (5.22h) n i+1/2/j/k = diag[max #( a i+1/2 )] A i+1/2 , 1 ^ (5.221) ni,j + i /2/k = diagtmax iHaj +1/2 )] Aj +1/2 , X* (5.22J) ni,j /k+1/2 = diag[max tf(a k+1/2 )] A k+1/2 ,

PAGE 88

77 Here &i +1/ j /k , &i,j+i,k> and ^i,j,k+l are Eq.(5.6) evaluated at (i+l,j,k), (i,j+l,k) and (i,j,k+l) respectively. Observe that Eq.(5.22) is the original Beam and Warming I i r algorithm [3] if 0i+i/2,j,k = "i,j+l/2,k = "i,j,k+l/2 = in Eq. (5.22h)-(5.22j) and $i+i/2 R i+l/2 in Eq(5.15) is replaced by the conventional fourt-order dissipation term. Therefore, the implemention of this conservative linearized ADI scheme into an existing central difference code (such as the one based on the Beam and Warming algorithm) is relatively simple. One has only to add the extra matrices fli+i/2,j,k' n f ni,j+l/2,k» anc * °i,j,k+l/2 for tne implicit operator and uses a more sophisticated dissipation term $i+i/2 R i+l/2 for tne explicit operator.

PAGE 89

CHAPTER VI THIN LAYER NAVIER-STOKES SOLVERS The basic governing equations and numerical algorithm used in this chapter have been adopted from the works of Steger [37], and Pulliam and Steger [36], The threedimensional Reynolds-averaged Navier-Stokes equations in generalized curvlinear coordinates are simplified by using the thin-layer assumption. These equations can be further reduced to a set of axisymmetric equations which govern azimuthal invariant flows. The resulting r; -invariant equations are solved by using an implicit, non-iterative approximate factorization finite difference scheme developed by Beam and Warming [3]. The axisymmtric thinlayer Navier-Stokes code, which was originally developed by the U.S. Army Ballistic Research Laboratory [26], [38], is modified to incorporate the TVD schemes derived in Chatpter V and to adopt a spatially variable time step procedure to enhance both the numerical accuracy and computational efficiency. In the following, a brief presentation is made to describe the assumptions and simplifications of the governing equations, the turbulent closure model adopted, and the implementation of the TVD schemes. 78

PAGE 90

79 6.1 Thin Layer Approximation In external high Reynolds number viscous flows, the effects of viscosity are usually concentrated near the no slip boundaries. Hence the flow profile near the solid surface must be accurately resolved, especially in the direction normal to the solid surface. In order to provide an adequate spatial resolution, grid lines near the solid boundary are exponentially clustered (so that at least one or two grid points will be within the viscous sublayer) . Without flow separation, it is usually assumed that only the viscous terms in the ($•-) direction normal to the solid surface are needed to be retained while the viscous terms in both the streamwise (£-) and spanwise (r,-) directions can be neglected. The resulting equations with this thin layer approximation can be written as follows : aq at dP a6 1 a § (6.1) + + + = dT d£ dr) 3f Re 9i where the viscous terms in the f -direction have been collected in the vector |, Detailed derivation for both the three dimensional full Navier-Stokes equations and the thin layer approximation in the generalized curvilinear coordinates are given in Appendix A.

PAGE 91

80 6 2 Azimuth al for n-) Invariant Equations Azimuthal-invariant equations are obtained from Eq. (6.1) by makinq use of followinq two restrictions: (l) all body qeometries are axisymmetric to the r? -direction and (2) the state variables and the contravariant velocities do not vary in the t? -direction. The equations resultinq from simplifications were derived by Nietubicz et al [26], [38] as shown in the followinq aq (6.2a) + 3E a6 l a§ dr 3£ a r Re d r where (6.2b) 1 — — __ __ P pU pu puU+£ x p q = J" 1 P v pW e _ • / E = j-i pvU pwU+e 2 p (e+p)U-£ t p pW pUW+f x p m 1 u f +m 2 f x 6 = J-1 pvW pwW+f 2 p ; 3 = = j-i m 1 w f +m 2 r z (e+p )W-f t pJ m 1 m 3 +m 2 (r x u+f A,

PAGE 92

81 (6.2c) ft = J-l and ° pV[R£(U-£ t )+R r (W-ft)] -pVR(V-i, t )-p/R *i = ^(r x 2 + fz 2 ). m 2 = (m/3) (r x u r + r z w f ), m 3 = (u 2 + v 2 + w 2 )j./2 + M Pr~ 1 (7-l)" 1 ( u = ft * ^xu + ^ z w, V = ft + ^yV/ w " ft fx u + r z w c K' 6.3 Turbulence Model The algebraic mixing length model developed by Baldwin and Lomax [39] is incorporated into the present thin layer formulation to model the effect of turbulence. This is a two-layer algebraic model. The inner layer is governed by the Prandtl mixing length with the Van Driest damping term; the outer layer follows Clauser's approximation. Computed vorticity is used in defining the reference mixing length required for the outer layer. The present turbulence model is appropriate for attached and mildly separated flow problems.

PAGE 93

82 6.4 Implicit Central Difference Algorithm The numerical scheme used here is an implicit, noniterative approximate factored, finite difference algorithm written in delta form, proposed by Beam and Warming [3]. The resulting finite difference equations can be expressed as n n n n n (6.3a) L^L r [ Aq ] = -At(S^t + S^ Re _1 « r § + ft ) D e , with the linear operators L^ and Lr as (6.3b) L£ [I + hS^ n Dil^], (6.3c) L f = [I + h5 r C n hRe-^J" 1 ft n J Dil r ]. Here h = At, the S's are typically three point secondorder accurate central difference operators and A and v are forward and backward-difference operators respectively. J is the Jacobian of the coordinate transformation. The Jacobian matrices A, C and ft result from the local time linearization and can be found in Appendix A. In the current study, the first-order accuracy in time and secondorder central difference in space are applied throughout the code Since central discretized factorization scheme was used, the high frequency error caused by the nonlinear terms may not be completely damped out, which can lead to

PAGE 94

83 numerical solutions with spurious oscillations. Therefore numerical dissipation terms denoted as Di and D e have been inserted into Eq. (6.3). The explicit fourth-order artificial dissipation term (6.4) D e = e e At J-i^A^V^) 2 + (A^Vj-) 2 ] J, is added to the right-hand side of Eq. (6.3) and the implicit second-order smoothing terms (6.5a) Di | ? = ei At J -1 (AV)^ J> (6.5b) Di | r = €± At J" 1 (AV) r J, are inserted into the respective implicit block operators. The implicit second-order smoothing terms are chosen to keep the left hand side factors of the block tridiagonal structure. The parameters e± and e e are chosen such that *i * 2e e . The solution algorithm now consists of two onedimensional sweeps, one in the fand the other in the $direction. Each step requires the solution of a linear system involving a block tridiagonal matrix, which is solved by a block LUD (lower-upper decomposition) solver. The resulting solution process is much more economical than the unfactored algorithm in computer storage and CPU time.

PAGE 95

84 Also steady state calculation can be achieved by integrating the governing equations from an arbitrary initial condition to a time asymptotic state. 6.5 Dissipation Model of TVD Scheme As mentioned in Chapter V, the major difference between the Beam-Warming scheme and the TVD scheme is that the TVD scheme requires more sophisticated dissipation terms. The dissipation terms of the Beam-warming scheme in Eq.(6.4) and Eq. (6.5) are replaced by the following expressions (6.6) D e = -H At(R i+1/2 $ i+1/2 Ri-i/ 2 *i-i/2 + R k+l/2*k+l/2 " R k-l/2*k-l/2)/ and (6.7a) Di \ ( = -h At (ni +1/2/k fli_ 1/2/k ), (6.7b) Di | r = -h At (ni /k+ i /2 ni,k-i/ 2 K where o 1+1 / 2 fk and n i,k+l/2 are defined similarly as in Eq. (5.22h) and Eq. (5,22i) except that the j index is ommitted.

PAGE 96

85 6.6 Grid System A schematic illustration of a Secant-Ogive-CylinderBoattail (SOCBT) projectile shape is shown in Figure 6.1. Specifically, the projectile under study is 3-caliber for secant-ogive part, 2 -caliber for the cylinder part and 1caliber 7 degree for boattail part which is further extended for another 1.77 caliber to meet straight sting. A grid generation code utilizing a hyperbolic differential equation [40] is employed to generate a network of constant £ and f lines in the physical x-z plane as illustrated in Figure 6.2. Corresponding uniform grids of £ and f in the computational space define a one to one mapping between points j,k in the physical plane to points j,k in the computational plane also shown in Figure 6.2. The grid spacing in the computational domain is chosen to be unity, that is A£ = l and Af =1. Once the projectile geometry is determined, the grid points along the body surface can be redistributed by using a cubic polynomial as the clustering function to make much denser grids near the high gradient region.

PAGE 97

86 a •p o
PAGE 98

87 Pysical Domain C = C 5=0 Computational Domain C = C / max V' //////// •9 5 J.k B 5 = Cmax 5=5 max /////// — C S=0 Figure 6.2 Mapping from physical domain to computational domain.

PAGE 99

88 For the projectile problem, since the outer boundary is unconstrained, it only needs to be sufficiently far away from the inner boundary. Such a grid generation problem is frequently encountered in external aerodynamic application. In the current study, the outer boundary is taken as 4 projectile lengths away from the projectile body. The " Ctype " grids are then generated by the hyperbolic formulation [41] . When solving the viscous flow problem it is necessary to have a very fine mesh normal to the projectile body near the wall, so that the velocity gradients within the boundary layer can be adequately resolved. The grid spacing can be smoothly controlled by using an exponential clustering function in the zdirection. The minimum grid spacing used near the wall in the code is 2xl0~ 5 diameters of projectile. 6.7 Numerical Boundary Conditions The general curvilinear coordinate transformations are made such that physical boundaries are mapped to boundaries in the computational domain. This makes the formulation and implementation of boundary conditions easier. For the axisymmetric projectile flow problem with an attached sting, unknown values of q on the boundaries are updated explicitly after each integration step. The numerical boundary conditions for the projectile calculations are as

PAGE 100

89 follows: (1) Along the outer boundary, line AB in Figure 6.2, free stream values are specified. (2) Along the line AD in Figure 6.2, the flow variables are obtained by a second-order extrapolation. Numerically this is done by using a three-point formula to approximate the first derivative with respect to £ . (3) Along the outflow boundary, line BC in Figure 6.2, all the flow variables are obtained by a two-point extrapolation. However, for the subsonic cases the pressure is specified along the outflow boundary(4) Along the inner boundary, line CD in Figure 6.2, the no-slip condition is employed for three velocity components (i.e., u,v, and w = 0) . The surface pressure is obtained from the normal momentun equation 3P/dn = 0, where n is the local normal to the solid surface. The surface density is obtained from an adiabatic wall assumption. Finally, total energy is computed from the equation of state. 6.8 Spatial-Varying Time Step For steady state flow problems, it is desirable to achieve the converged solution as quickly as possibble. The convergence to a steady-state solution can be accelerated by using spatially varying time steps [4]. For a fixed time step, the Courant number is not uniform since the grid spacings vary from very fine to very coarse in the

PAGE 101

90 whole flow field. The use of a spatial varying time step can thus be interpreted as an attempt to use a more uniform Courant number throughout the field. The choice of this spatial varying time step is governed by the formula At ref (6.8) At = 1 1 + (J)^ where the J is the Jacobian of transformation.

PAGE 102

CHAPTER VII RESULTS AND DISCUSSION 7.1 Numerical Implementation One objective of the present study is to compare the conventional Beam-Warming scheme with the TVD scheme using different flux limiters in terms of accuracy, stability and efficiency. In this chapter the turbulent transonic flow field about a secant-ogive-cylinder-boattailed (SOCBT) projectile model is adopted as the test problem. An axisymmetric thin-layer Navier-Stokes code developed at the U.S. Army Ballistic Research Laboratory [26,38] has been modified in the present study. One major difference between the original and modified codes is the treatment of the numerical dissipation terms. This approach of implementing the TVD scheme is suggested by Yee [42], who interpreted the TVD algorithms as a form of nonlinear numerical dissipation. This viewpoint can greatly reduce the work of incorporating the TVD schemes into the existing computer code. In the original version of the code, as reported in [26], the BeamWarming scheme was adopted. Based on this framework, the TVD schemes can be incorporated into the original computer code by adding a more sophisticated numerical dissipation model. With 91

PAGE 103

92 this improved numerical dissipation mechanism, one can strive both to eliminate the non-physical oscillations commonly encountered in the numerical solution yielded by the BeamWarming scheme and to satisfactorily resolve the sharp gradients present in the flow field. In the present investigation, three variants of the flux limiters have been assessed in detail; two of them are designed for the upwind TVD scheme, and the third one is for the symmetric TVD scheme. Specifically, the flux limiter, Eg(5.19b), of the upwind TVD scheme suggested by van Leer is called 'Limitl'; the one which combines Roe's superbee limiter, Eq. (5.19c), (in the linear field) and van Leer's limiter ,Eq. (5.19b), (in the nonlinear field) is called 'Limit2', The flux limiter of the symmetric TVD scheme, Eq. (5.21c), is called 'Limit3'. Both the original and modified thin-layer Navier-Stokes codes have been used to compute transonic turbulent flow over the SOCBT projectile model at zero angle of attack. The test cases considered here are the SOCBT projectile with (a) Me = 0.96, (b) Mo, = 0.91, and (c) M*, = 1.2. Here M*, is the free stream Mach number. Figure 7.1 shows a 90x60 " C-type" grid for a flow domain of 4 projectile lengths with special clustering on both the ogive-cylinder and cylinder-boattail junctures. Moreover, an exponential clustering function is employed such that the minimum grid spacing chosen above the

PAGE 104

93 Figure 7.1 The 90x60 "C-type" grid for the SOCBT projectile with sting. .' 'I' 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 [it 1 1 1 1 1 1 1 1 1 1 1 1 Beam-Warming scheme 90x60 "C-type" grid I.M9MJ S.^m-H LSXt-IH (.MM lU-« Figure 7.2a The surface pressure coefficients of Beam-Warming scheme for the SOCBT projectile M„ = 0.96, and Re = 7.7xl0 5 .

PAGE 105

94 projectile surface in the f -direction is 2xl0~ 5 diameters of the projectile in order to adequately resolve the viscous sublayer. For the original thin-layer code, the parameters of user-controlled smoothing terms are set to e^ = 2c e = 4At. For the modified code, the coefficient of e in Eq.(3.11) is set to zero; also the spatially varying time step, Eq. (6.7), can be used to further enhance the convergence rate. No special effort was devoted to vectorizing the program for either version of the code. Hence the CPU time presented for the following cases should be viewed only as reference. All the numerical computations were run on the CRAY-XMP/48 supercomputer at the National Center for Supercomputing Applications (NCSA) , at the University of Illinois at UrbanaChampaign. 7.2 Results of Flow of H = 0.96 Figure 7.2a shows the converged solution for the surface pressure coefficient C p calculated by the original BeamWarming scheme. The symbols " o " on the plot indicate the experimental data obtained by Kayser et al [27], while the solid line represents the numerical solution. As one can see in Figure 7.2a, the numerical result compares favorably with the experimental data except for the peak point located in the shock region above the middle of cylinder segment. The corresponding comparision for the results obtained

PAGE 106

95 with the TVD schemes are presented in Figure 7.2(b,c,d). Overall, the surface pressure coefficient yielded by both the Beam-Warming and TVD schemes are of similar quality on the present 90x60 grid system. A comparison of the pressure and shear drag coefficients predicated by the four schemes is made in Table 7.1. The nondimensionalized pressure drag coefficients, defined as Dp = (pressure force) /(p M a*, 2 L 2 ) obtained by all four schemes are very close, with less than 1 percent difference between any two predictions. The nondimensionalized shear drag coefficient, defined as D f = (shear force)/ (0.5 Pa> V^ 2 ) shows more discrepancies among the predictions. The BeamWarming scheme predicts the smallest value of Df, while the symmetric TVD scheme predicts the largest shear drag among all the schemes under study. The difference between the largest and the smallest predicted values of Df is more than two percent. Although there are no experimental measurements available for a definite assessment with respect to the drag coefficients, the results shown in Table 7.1 broadly confirm the effects of the TVD schemes on the numerical solution.

PAGE 107

z UJ u t Ui o g at Q. X/D Figure 7.2b The surface pressure coefficients of upwind TVD scheme with Limit 1 for the projectile with M*, = 0.96, and Re = 7.7xl0 5 using a 90x60 C grid. z UJ u t UJ O u «n UJ tc a. X/0 Figure 7.2c The surface pressure coefficients of upwind TVD scheme with Limit2 for the projectile with Ma, = 0.96, and Re = 7.7xl0 5 using a 90x60 C grid.

PAGE 108

97 o t i « ! MflCH NO. = 0.96 o exper inonial data nunerJcal volue X/0 Figure 7 . 2d The surface pressure coefficients of symmetric TVD scheme with Limit3 for the projectile with Me = 0.96, and Re = 7.7x10 s using a 90x60 C grid.

PAGE 109

98 In terms of the computational efficiency, Figure 7.3 shows the convergence history of the four schemes under study. The symmetric TVD scheme has the fastest convergence rate. Furthermore, while there exist differences in the rate of convergence among the three TVD schemes, Figure 7.3 clearly indicates that their perforamces are all much superior to that of the Beam-Warming scheme. The other advantage of the TVD schemes is the fact that the user does not have to resort to empirical rules to choose the damping parameters. In the Beam-Warming scheme, the proper choices of the free parameters such as q and e e are not always clear; a trial-and-error approach is usually needed. In this regard, the TVD schemes can be viewed as a more rigorous and sophisticated approach in determing the desirable numerical dissipation. Table 7.2 presents more information of the computational aspects of the four schemes. Only fixed time step is used by the Beam-Warming scheme, whereas TVD schemes are used both with and without the provision of the spatially varying time step. On the basis of the CPU time per iteration, the BeamWarming scheme is the most efficient. This is clearly expected since the TVD schemes need to make a substantial amount of extra computations to eliminate the spurious oscillations and to produce high profile resolutions. Among the TVD schemes, the symmetric TVD scheme has the best performance in terms of CPU time per iteration. Combined with

PAGE 110

99 i i i — i — |— i — i i i i i | i i i r o.wtt-n tr 2: cc o z CM anc-B 0. JOIt-B MflCH NO. = 0.96 I I ' ' I I I L2MiO% O.fMfrfCM ITERATIONS V Figure 7.3 Convergence history of different schemes (a) Beam-Warming scheme (b) Upwind TVD scheme with Limit2 (c) Upwind TVD scheme with Limit 1 (d) Symmetrtic TVD scheme with Limit3

PAGE 111

100 the superior convergence rate, the symmetric TVD scheme is the most efficient scheme overall. But nevertheless, Table 7.2 demonstrates that with the identical grid distribution, all of the TVD schemes can yield the steady state solution with much less CPU time than the Beam-Warming scheme. With the aid of the spatially variable time steps, significant improvements in convergence rate have been obtained. The symmetric TVD scheme is the most attractive one for the present turbulent transonic flow calculations. Table 7.1 Pressure and Shear Drags for Different Schemes of 90x60 Grid at *!«, = 0.96 Scheme Limiter Pressure Drag x 10~ 3 Shear Drag x 10" 3 B-W 40.96 27.88 UPTVD Limit 1 40.56 28.22 UPTVD Limit2 40.67 28.14 SYMYVD Limit3 40.72 28.52

PAGE 112

101 Table 7.2a L 2 -Norm Residual and CPU Time for Different Schemes with 90x60 Grid at Mc = 0.96 Using Fixed Time Step Scheme Limiter L2 Norm Residual NO. Of Iter. CPU per Iter. Total CPU (min) B-W 1 9.6xl0~ 5 6000 0.73 72.14 UPTVD 2 Limit 1 5.8X10" 5 1700 0.92 25.00 UPTVD Limit2 7.1X10" 5 1700 0.96 27.00 SYMTVD 3 Limit3 7.2X10" 5 1300 0.85 17.65 1 : Beam-Warming Scheme 2 : Upwind TVD Scheme 3 : Symmetric TVD Scheme Table 7.2b L 2 -Norm Residual and CPU Time for Different Schemes with 90x60 Grid at M, = 0.96 Using Spatially Varying Time Step Scheme Limiter L2~Norm Residual No. of Iter. CPU per Iter. Total CPU (min) UPTVD Limit 1 6.3xl0~ 5 1000 0.90 15.18 UPTVD Limit2 5.1X10" 5 1400 0.95 22.35 SYMTVD Limit3 7.2X10"" 5 700 0.81 10.72 7.3 Grid Density and Numerical Accuracy In Figure 7.2, satisfactory and comparable numerical accuracies among the four schemes were obtained on a 90x60

PAGE 113

102 grid system. In order to assess the effectiveness of the TVD schemes and the relative merits between them and Beam-Warming scheme, a series of computations has been conducted by using a progressively coarser grid distribution. First, the number of grid points along the streamwise, i.e., £-, direction is reduced from 90 to 60. That is, the grid system now employs 60x60 nodes instead of 90x60 nodes. Figure 7.4 presents the comparisons between the experimental results and numerical ones in the same flow condition on this grid system among the four schemes. Compared to Figure 7.2, the solutions obtained by the TVD schemes on both the 90x60 and 60x60 grid system are essentially undistinguishable in term of the pressure coefficient. The Beam-Warming scheme, on the other hand, shows degradation in accuracy on the coarser grid system, notably in the shock region. Next, the number of grid points along the transverse, i.e., f-, direction is reduced from 60 to 40. Figure 7.5 compares the prediction and measurements of the pressure coefficient on a 90x40 grid system for the Beam-Warming scheme and symmetric TVD scheme. Overall, the Beam-Warming scheme fails to capture the shock adequately and its accuracy is much inferior to the results obtained by using the TVD schemes. Again all three variants of the TVD scheme perform comparably well, although only one of them is presented in Figure 7.5. The impact of the reduction of grid

PAGE 114

103 UJ G o MfiCH NO. = 0.96 o ixpcrinenlal doia numerical value x/o Figure 7.4a The surface pressure coefficients of Beam-Warming scheme for projectile with M<*> 0.96, and Re = 7.7x10 s using a 60x60 C grid. The L2~norm residual of 5.16xl0~ 5 can be reached in around 5000 time steps. i i i — MfiCH NO. = 0.96 o experinerrlal daic numerical value X/0 Figure 7.4b The surface pressure coefficients of upwind TVD scheme with Limitl for projectile with M» = 0.96, and Re = 7.7xl0 5 using a 60x60 C grid. -5 The L 2 -norm residual of 3.9xlO~ D can be reached in around 1500 time Steps.

PAGE 115

104 MflCH NO. = 0.96 o experimental data nurerJcal value Figure 7.4c x/o The surface pressure coefficients of upwind TVD scheme with Limit2 for projectile with M» = 0.96, and Re = 7.7x10 s using a 60x60 C grid. The Lo-norm residual of 8.5x10-° can be reached in around 1900 time steps, UJ i — -i 1 1 — HfiCH NO. = 0.96 o experimental data nurerJcal value X/D Figure 7.4d The surface pressure coefficients of symmetric TVD scheme with Limit3 for projectile with M» = 0.96, and Re = 7.7xl0 5 using a 60x60 C grid. The L 2 -norm residual of 3.4xl0 -5 can be reached in around 1100 time steps.

PAGE 116

t o g 105 1 1 — MflCH NO. = 0.96 o experinerrtal data nurerical value a.mm-m x/o Figure 7 . 5a The surface pressure coefficients of Beam-Warming scheme for projectile with M« = 0.96, and Re = 7.7xl0 5 using a 90x40 C grid. The L2-norm residual of 2.59xl0" 6 can be reached in around 3000 time steps. Q MflCH NO. = 0.36 o ex peri net-rial data nurerical value Figure 7 . 5b x/o The surface pressure coefficients of symmetric TVD scheme with Limit3 for projectile with M« = 0.96. and Re = 7.7xl0 5 using a 90x40 C grid. The L 2 -norm residual of 4.8xl0 -r> can be reached in around 1100 time steps.

PAGE 117

106 size on the numerical accuracy at the present level of grid size is not visible for the TVD schemes. Further reductions of the number of grid points have also been performed. Figures 7 . 6 and 7 . 7 show the comparisons of numerical results and measured data of the three TVD schemes on the 80x40 and 70x40 grid system respectively. Again, the three schemes all appear to be able to maintain good accuracy in both cases. In summary, by reducing the grid size from 90x60 to 70x40, the TVD schemes are shown to be both robust and accurate enough to yield a steady-state solution. Table 7 . 3 presents some comparisons of the performance among the TVD schemes on the 70x40 grid in the computational aspect. Overall the symmetric TVD scheme is fairly efficient in term of iterations and CPU time per iteration. However the convergence rate for symmetric TVD scheme has slightly degraded as shown in Figure 7.8. Table 7 . 4 compares the predicted pressure and shear drags on the 70x40 grid system. It appears that results yielded by upwind TVD with Limitl and the symmetric TVD with Limit3 closely follow each other, while that by the upwind TVD scheme with Limit2 predicts smaller drag forces.

PAGE 118

z fc UJ o u I
PAGE 119

108 o MRCH NO. = 0.96 o ex peri rental doia nunerJcal value X/0 Figure 7 . 6c The surface pressure coefficients of symmetric TVD scheme with Limit3 for projectile with M<*> = 0.96, and Re = 7.7xl0 5 using a 80x40 C grid. The L2~norm residual of 5.4xl0~ 5 can be reached in around 1100 time steps.

PAGE 120

109 z a 4T> MflCH NO. = 0.96 o experinenial data nunerjcal value X/0 Figure 7.7a The surface pressure coefficient of upwind TVD scheme with Limitl for projectile with M^ 0.96, and Re = 7.7xl0 5 using a 70x40 C grid. 8 X/0 Figure 7.7b The surface pressure coefficient of upwind TVD scheme with Limit2 for Projectile with Me = 0.96, and Re = 7.7xl0 5 using a 70x40 C grid.

PAGE 121

\ — I t 8 no MflCH NO. = 0.96 o experimental daia nunerical value X/0 Figure 7.7c The surface pressure coefficient of symmetric TVD scheme with Limit3 for projectile with Mo — 0.96, and Re = 7.7xl0 5 using a 70x40 C grid.

PAGE 122

Ill o. x»c-«i -1 — I — I — I I I — I I I — I — I I I I I II — I I I — I— I — I I MRCH NO. = 0.96 O.KK-B .UCEfOi _i I L J I I I I.IKCiO* D.I ITERATIONS Figure 7.8 The convergence history of different TVD schemes at M* = 0.96, using a 70x40 grid. (a) Upwind TVD scheme with Limit2 (b) Upwind TVD scheme with Limitl (c) Symmetric TVD scheme with Limit 3

PAGE 123

112 Table 7.3 L2-Norm Residuals and CPU Time for Different Schemes of 70x40 Grid at M*, = 0.96 Scheme Limiter L2~Norm Residual No. of Iter. CPU per Iter. Total CPU (min. ) UPTVFD Limitl 4.2X10" 6 1500 0.44 10.90 UPTVD Limit2 1.4X10" 5 1900 0.49 17.27 SYMYVD Limit3 3.2X10" 5 1100 0.42 7.83 ! Table 7.4 Pressure and Shear Drags for Different Schemes for 70x40 Grid at M*, = 0.96. Scheme Limiter Pressure Drag xl0~ 3 Shear Drag xl0~ 3 UPTVD Limitl 43.42 28.09 UPTVD Limit2 41.57 27.68 SYMTVD Limit3 43.20 28.05 7.4 Results of Flow of M^ = 0.91 In this case the transonic speed is reduced from M^ = 0.96 to Me = 0.91. All the calculations were conducted on the 70x40 grid system. The computational aspects of the TVD schemes are described in Table 7.5. The computed pressure and shear drags are shown in Table 7.6. Figures 7.9(a,b,c) show that the agreements of computed surface pressures with

PAGE 124

113 UJ O i. i MRCH NO. a 0.91 o experimental data numerical value X/D Figure 7 . 9a The surface pressure coefficients of upwind TVD scheme with Limit 1 for projectile with M,, = 0.91, and Re = 7.7xl0 5 using a 70x40 C grid. t O o O in MRCH NO. = 0.91 o experimental data numerical value X/D Figure 7.9b The surface pressure coefficient of upwind TVD scheme with Limit2 for projectile with M» = 0.91, and Re = 7.7xl0 5 using a 70x40 C grid.

PAGE 125

114 u fc UJ o 1 UJ ^ 1 1 — HfiCH NO. = 0.91 o experinenial doio numerical value X/0 Figure 7.9c The surface pressure coefficients of symmetric TVD scheme with Limit3 for projectile with Me = 0.91, and Re = 7.7xl0 5 using a 70x40 C grid.

PAGE 126

115 experimental data are fairly good. Figure 7.10 shows the convergence history of the TVD schemes with three different flux limiters. The convergence rate is degraded when the upwind TVD scheme with Limit2 is used. Roe's " superbee " limiter used in the linear field tends to produce more compressive profiles than other limiters. The results obtained here suggested that the convergence rate will degrade when too compressive limiter is used for the flow with weak shocks. Table 7.5 I^-Norm Residuals and CPU Time for Different Schemes of 70x40 Grid at M^, = 0.91 Schemes Limiter L2 -Norm No. of CPU per Total CPU Residual Iter. Iter. ( Min.) UPTVD Limitl 1.3xl0~ 5 1900 0.47 14.75 UPTVD Limit2 3.4xl0~ 4 1300 0.49 12.29 SYMYVD Limit3 5.1X10" 6 1500 0.42 10.61 Table 7.6 Pressure and Shear Drags for Different Schemes of 70x40 Grid at M*, = 0.91. Schemes Limiter Pressure Drag x 10~ 3 Shear Drag x 10" 3 UPTVD Limitl 21.82 25.25 UPTVD Limit2 18.80 23.07 SYMTVD Limit3 20.97 25.29

PAGE 127

116 BJ O.MK-« a — i cr> £ £ o.Mtc-n o CM -J 6.MC-M I I I I — I — 1—1 — I — I — I — I — 1—1 — I I I I I I I MflCH NO. = 0.91 0. 10K-B ' — ' 1 — I 1 — I — ' — ' — I — I — I — I — — i 1 — I — 1 — i I I i | I I Ltn«X 0.KCC.O* I.ZMtO* O.tMMI «.«S«iO» D.SKK01 ITER ATI ONS Figure 7.10 The convergence history of different TVD schemes at M^, 0.91 using a 70x4 grid. (a) Upwind TVD scheme with Limit2 . (b) Upwind TVD scheme with Limitl. (c) Symmetric TVD scheme with Limit3 .

PAGE 128

117 7.5 Result of Flow of 1L = 1.2 In this case, as the Mach number is increased from M,*, = 0.96 to Ma, = 1-2, computations on the 70x40 grid system have been conducted for the three TVD schemes. Figures 7.1l(a,b,c) show that surface pressure coefficients computed by using the TVD schemes are in excellent agreement with the measured data. All three TVD schemes appear to produce almost identical results. Figure 7.12 shows the convergence histories of the three TVD schemes. The convergence rate of the symmetric TVD scheme has shown difficulties in reducing the residual to a very low level whereas the other two schemes can rapidly reduce the residuals throughout the whole calculation. It is noted that for all the calculation with TVD schemes, the coefficient of entropy correction e is assigned to be zero in the current code. Yee has pointed out that the smaller the € being used, the slower is the convergence rate of the symmetric TVD scheme [42]. Although the convergence rate could be further improved if a positive value of e is used, it was decided not to pursue this in the present work since one of the major attractions of the TVD scheme is to be free from arbitrarily choosing the adjustable parameter values. The pressure drag and shear drag of this case are shown in Table 7.7. The predictions of the two drag coefficients are very close among the three schemes.

PAGE 129

118 z *•* o " g UJ B 1 I | i n MRCH NO. = 1.2 o experimental doia nunerJcal value -Q.CQj£*m " LAAAiLMtlAA.. . . . . 1 I.MDKQ O.MOffOI t.j o tawi «; L x/o Figure 7 . 11a The surface pressure coefficients of upwind TVD scheme with Limitl for projectile with Ma = 1.2, and Re 7.7xl0 5 using a 70x40 C grid. t 11 » ' * I ' I ' | ' « *^'l UH | ltl1 TT ' ' i i n i -r-r i i i r r ^ MRCH NO. = 1.2 o experjnenla] doia nunerJcal value aa
PAGE 130

u ! «n A 119 't 1 * wm MflCH NO. = 1.2 o experimental data nurerJcal value X/D l«« O.nWTKX IM*. raiHOI *.M<» 0.TSU1 Figure 7.11c The surface pressure coefficients of symmetric TVD scheme with Limit3 for projectile with Me, = 1.2, and Re = 7.7xl0 5 using a 70x40 C grid.

PAGE 131

120 a z § 1 — i — " — i — i — i — i — i — i — i — r -i — I — i rMRCH NO. = 1.2 O.MC-W O.MC-M 0.WK-B : -i — i i i_ O.IOOCtOt •.IBCiOt -I 1 l__l i I MCiOt D.IOCfMM ITER ATI ONS Figure 7.12 The convergence history of different TVD schemes at M. = 1.2 using a 70x40 grid (a) upwind TVD scheme with Limit2, (b) upwind TVD scheme with Limitl and (c) symmetric TVD scheme with Limit3.

PAGE 132

121 Table 7.7 Pressure and Shear Drags for Different Schemes of 70x40 Grid at Ma, = 1.2. Scheme Limiter Pressure Drag xlO"* 3 Shear Drag xlO" -3 UPTVD Limitl 99.51 43.61 UPTVD Limit2 99.50 43.52 SYMTVD Limit3 96.66 43.51 The results presented in the previous section clearly illustrate that both symmetric and upwind TVD schemes are able to capture shocks accurately and efficiently. The numerical solutions obtained at different transonic speeds show that as the shock strength increases, the TVD schemes appear to produce more accurate results.

PAGE 133

CHAPTER VIII SUMMARY AND CONCLUDING REMARKS In the present dissertation, extended derivation of the TVD schemes in three-dimensional generalized curvilinear coordinate systems has been presented. The derived schemes are then simplified and implemented into a time-dependent axisymmetric thin-layer Navier-Stokes flow code by improving the treatment of the numerical dissipation terms to make the resulting overall scheme consistent with the TVD formulas. A spatially varying time stepping procedure has also been incorporated into the code to accelerate the convergence rate of the steady state flow computations. This improved version of the code was used to predict transonic turbulent flows of Ma, = 0.91, 0.96 and 1.2 over an axisymmetric SOCBT projectile at zero angle of attack. The numerical results show good agreement with available experimental surface pressure coefficient data. The use of the TVD schemes requires more CPU time per time step than the conventional Beam-Warming scheme. However, substantial improvements of the convergence rate have been observed with all three variants of TVD scheme. Hence, overall the TVD schemes are more economical to apply for the flows studied here. Furthermore, equally accurate 122

PAGE 134

123 solutions can be obtained with the TVD schemes on less number of grid points compared to those with the Beam-Warming scheme. In terms of the relative merits among the three variants of the TVD scheme, comparable accuracies are yielded by both the upwind and symmetric schemes. However, it is found that for the flow with a weak shock (K^ = 0.91), the flux limiter proposed by Roe for the upwind TVD scheme can experience some difficulties in convergence. On the other hand, for the flow with stronger shock (M*, = 1.2), the symmetric TVD scheme can experience difficulties in the later stage of convergence. It is concluded that although the flux limiter proposed by van leer for the upwind TVD scheme may not be always the most efficient one among the three variants studied here, its convergence appears less affected by the incoming flow Mach number. The robustness of this scheme makes it a favorable choice in practical flow computation.

PAGE 135

APPENDIX A NAVIER-STOKES EQUATIONS AND THIN LAYER APPROXIMATION The strong conservation law form of the three-dimensional Navier-Stokes equations in Cartesian coordinates can be written in the flux vector form as (A-l) dq 3E 3F 3G dt ax 3y dz 3 Ey a Fy + + ax ay aG v az where (A-2) q = p P u pv pW e E = pU p\JL 2 +p p\XV puvr u(e+p) , F = P v pVU pv 2 +p pVW v(e+p) , G = pv pWV pw 2 +p w(e+p) E v = Re -1 T xx T yx T zx Px , F v = Re -1 r xy ryy T zy 3 _ Gv = Re -1 T— —1 r xz r yz r zz Pz 124

PAGE 136

with 125 (A-3) where Tyy T zz T xy r xz = r yz = £x = £y = = r A(U X + V y + W Z ) + 2/iU x , A(U X + Vy + W Z ) + 2/iVy, A(U X + Vy + W 2 ) + 2AiW z , yx = /* (Uy + v x ) , T ZX = M(U Z + W x ) , r zy = A«(v z + Wy) , Ur XX + Vr xy + Wr xz + / iPr1 (7-l)" 1 a x C 2 , UTy X + VTyy + WT y z + ft Pr" 1 ( 7 1 ) ~ 1 d yC 2 , Ur zx + Vr zy + wr zz + ^Pr -1 ( 7 -l) _1 a Z C 2 , p is the density, made dimensionless by p m , c is the local speed of sound and c 2 = 7p/p , u, v, w are Cartesian velocity components nondimensionalized by speed of sound c, e is the total energy nondimensionalized by p m c 2 , p is the pressure ( = (7-1) [e 0.5p (u 2 +v 2 +w 2 ) ] ) , 7 is the ratio of specific heats, A* is the dynamics viscosity, A is from Stokes's hypothesis, = -2/3/*, k is the coefficient of thermal conductivity, Re is the Reynolds number, and Pr is the Prandtl number.

PAGE 137

126 Transformed Navier-Stokes Equations To enhance numerical accuracy and efficiency and to handle boundary conditions more easily, the Navier-Stokes equations can be transformed from Cartesian coordinates to general curvilinear coordinates where (A-4) r = t ( = l(x, y, z, t) v = r, (x, y, z, t) T = T (x, y, z, t) The resulting transformed equations can be written in nondimensional form as (A-5) 8r q + d£ d 8 3 . (E A E V ) + (F F v ) + (6 6 V ) dr, 0, ar where q = J-l (A-6) pU pW e E = J-l P U pUU + £ x p pvU + £ y p pWU + £ Z P (e+p)U £ t P

PAGE 138

I = J-l 127 pV puV + ry x p pW + r/yp pWV + »? z p (e+p)V r/ t P 6 = j-i pW pUW + f x p pVW + f y p pWW + f z p (e+p)W rtP and (A-7) u = it + C X u + CyV + £ z w, v = »?t + »?x u + ^yV + »7z w » W = ft + Tx u fyv + Tz w where U, V, W are contravariant velocity components. The viscous flux terms are given by £ v = (JRe) -1 fx r xx + £y r xy + ?z r xz £x r yx + fy r yy + £z r yz X T ZX + c yrzy + i Z T ZZ ZxPx + £y0y + £ z £z (A-8) £ v = (JRe)" 1 ^x T xx + , ?y r xy + '7z r xz »?x r yx + r iy T YY + r iz T yz lx r zx + »?y r zy + r lz T zz ixP: TlyPy + T)zPz &v = (JRe) -1 fx T xx + fy r xy + fz T xz fx r yx + fy T yy + fz r yz ^x r zx + Ty^zy + fz T zz Tx^x + fy/Sy + f 2 £ z

PAGE 139

128 where the components of the stress tensor and heat-flux vector in nondimensional form are given in Eq. (A-3) . The Cartesian derivates of the stress tensor can be expanded in £ , r? , and f spaces via the chain rule, for example u x = £x u £ + Vx^r, + fx u f • The metric terms are defined as £x = J(y^z r y r z^) r, x = J(z ? y r y^z r ), £ y = J(z,,x r x^z r ) n y = J(x^z r x z z£), £z = J(x^y r y^x^) r, z = J(y^x r X|y r ), (A-9) r X = J(Y?Zr, Z£Y„) ? t = "X r e X ~ Yr^y ? Z r £ z , Ty = J(X^Z^ X^Z^) vt = -Xtijx yrrjy z r r ?2 , and ^ = J( ^** " ***** rt = ~ x ' r * " Yrfy z r r z , = X£ Yr/ z r + x r y ?2r? + x^y r z^ x^z,, x^y^z r x r y„ Z £ J-l = Thin-Laver Approximation In high-Reynolds-number flows, the viscous effects are confined to a thin layer near no-slip boundaries. In most cases, there are only enough grid points to resolve the gradients normal to the body by clustering the grid in the normal direction. As a result, the viscous derivatives associated with the £ and r, directions are dropped and the

PAGE 140

129 terms in the f -direction are retained. With this assumption, Eg. (A-5) simplifies to (A-10) a T q +a^E +a^F +a r d = Re _1 a r §, where (A-ll) 3 = J -1 A»m!U r + ( M /3)m 2 f x /im^ + ( M /3)m 2 r y /imiw r + (M/3)m 2 r z / im 1 m 3 +(M/3)m 2 (fx u+ ry v+ fz w ) where (A-13) m l = fx + fy + fz/ m 2 = fx u r + **y v r + ^z w f. m 3 = (u 2 +v 2 + w 2 )/2 + Pr _1 ( 7 -1) (c 2 ) f

PAGE 141

APPENDIX B BEAM-WARMING IMPLICIT BLOCK ADI ALGORITHM AND ARTIFICIAL DISSIPATION The finite-difference algorithm due to Beam and Warming applied to Eq. (A-10) results in the following approximate factorization : (B-l) L e L,, L r Aq = (RHS) n , where L^ = (I + hs^k n + DjJ^), L„ = (I + hS„B n + Di|„) f L r = (I + h* r C n hRe1 * r J1 M n J + Dil r ), (RHS) n = -At(^E n + ^F n + 5 r d n -Re-ifij-finj D e , and where 6 is the central-difference operator and A and v are forward and backward-difference operators, and h = At corresponds to the first-order time-accurate Euler implicit scheme. DjJ^, Di | n and Di | jare the second-order implicit and D e is the fourth-order explicit smoothing operator. The Jacobian matrices A n , B n , C n and M n are obtained by 130

PAGE 142

131 linearizing the flux vectors £ n , $ n & n and il n in time such that (B-2) fcn+1 = £n + £n A £ + (h 2 ), f.n+1 . fn + fcn A £ + (h 2 ) , £n+l = £n + £n A * + (h 2 ), Re -l§n+l = Re -l[gn + j-lftn^nj + (h 2 ) , where k = dt/dq, 6 = d$/dq, £ = dt/dq, and ft = 3§/dq are the flux Jacobians. The Jacobian matrices are i or I or 6 = (B-3) «y i K X 4> 2 -U.9

k z $ 2 -w9

2 )e here «x(7ep 1 ~ 2 ) -( 7 -l)u0 Kyii/ep' 1 -* 2 ) I KzC-yep" 1 -^ 2 ) -( 7 -l)v0 | -( 7 -l)w* = K X U + Ky + /C Z W, .2 = = ( 7 -l) (u^ + v^ + w^)/2,

PAGE 143

132 with k = £, or i) or f for A, §, or, 6 respectively, The viscous flux Jacobian is (B-4) ft = j m 2 i "l^f (p -1 ) <*2 3 f (P _1 ) «3 a f (P _1 ) m 31 a 23 r (p1 ) Q 4 a r (p _1 ) a 5 a r (p -1 ) m 41 "3 a f (p _1 ) «5 a r (p -1 ) «6 a r (p-i) m 51 m 52 m 53 m 54 "0 a {(p _1 ) where m 21 = ttl a r (u/p)-a 2 a f (v/p)-a 3 d f (w/p) , m 31 = -<*2 a f (u/p)-a 4 3 r (v/p)-a 5 3 r (w/p) , m 41 = -«3*J (u/p)-a 5 3 r (v/p)-a 6 a f (w/p ) , m 44 = a 4 3 f (p -1 ) , m 51 = "O^f [-(e/p 2 )+(u 2 +v 2 +w 2 )/p] ai a r (u 2 /p)-a 4 a r (v 2 /p)-a 6 a r (w 2 /p) -2a 2 a f (uv/p)-2a 3 a r (uw/p)-2a 5 a r (vw/p ) , m 52 -a a f ( u /p) " ^21/ m 53 = -«oa$-(v/p) m 31 , m 54 " _ «O a f( w /p) " »41/ a = 7/aPr l (f fy 2 + fz 2 ), 2 _ r 2 _ «1 = M[(V3)r x Z fy" ~ fz z ]/ Q 2 = (^/3)f x fy>

PAGE 144

133 <*3 = (^/ 3 )rxfz# " y •4 = ^tTx 2 + (V3)fy 2 + f z 2 ], <»5 = (M/3)r x rzf «6 = M[f x 2 fy 2 + (V3)r z 2 ]. Artificial Dissipation The numerical solution of the Euler/Navier-Stokes equations for flows involving captured shock waves require the use of numerical or artificial dissipation. The reason for adding artifical dissipation is to control strong nonlinear effects such as shocks. The block algorithm which was given in Eq. (B-l) uses a constant coefficient dissipation model. In this approach, an explicit fourth-order dissipation (B-5) D e q n = £e At ^[(V^) 2 + (V^A,,) 2 + (V r A r ) 2 ] Jqn, is added to the right-hand side of the equation and an implicit second-order dissipation (B-6) Di| k = ei j-l A k v k J, * f ,*,t. is inserted into the respective implicit block operators.

PAGE 145

REFERENCES [1] P.D. Lax and B. Wendroff, "Systems of Conservation Laws," Comm. Pure Appl. Math., Vol. 13, 1960, pp. 217-237. [2] R.W. MacCormack, "The Effect of viscosity in Hypervelocity Impact Cratering," AIAA Paper 69-354, Cincinnati, Ohio, 1969. [3] R. Beam and R.F. Warming, "An Implicit Factored Scheme for the Compressible Navier-Stokes Equations," AIAA Paper 77-645, 1977. also AIAA J. Vol. 16, 1978, pp. 393-401. [4] T.H. Pulliam and J. Steger, "Recent Improvements in Efficiency, Accuracy and Convergence for Implicit Approximate Factorization Algorithms," AIAA Paper 85-0360, 1985. [5] A. Harten, "A High Resolution Scheme for the Computation of Weak Solutions of Hyperbolic Conservation Laws," J. Comp. Phys., Vol. 49, 1983. pp. 357-393. [6] A. Harten, J.M. Hyman and P.D. Lax, "On Finite-Difference Approximations and Entropy Conditions for Shocks," Comm. Pure Appl. Math., Vol. 29, 1976, pp. 297-322. [7] R. Courant, E. Isaacson and M. Rees, "On the Solution of Nonlinear Hyperbolic Diffential Equations by Finite Differences," Comm. Pure Appl. Math., Vol. 5, 1952, pp. 297322. [8] S.K. Godunov, "Finite Difference Method for the Numerical Computation of Discontinuous Solutions for the Equations of Fluid Dynamics," Math. Sbornik, Vol. 47, 1959, pp. 271-306. Also Cornell Aeronautical Lab. translation. [9] J. P. Boris and D.L. Books, "Flux-Correct Transport I SHASTA : A Fluid Transport Algorithm That Works," J. Comp. Phys., Vol. 11, 1973, pp. 38-69. [10] S.T. Zalesak, "Fully Multidimensional Flux Correct Transport Algorithms for Fluid," J. Comp. Phys., Vol. 31, 1979, pp. 335-362. [11] A. Harten, "On a Class of High Resolution TotalVariation-Stable Finite Difference Schemes," SIAM J. Numer. Anal. Vol. 21, 1984, pp. 1-23. 134

PAGE 146

135 [12] H.C. Yee, R.F. Warming and A. Harten, "Implicit Total Variation Diminishing (TVD) Schemes for Steady-State Calculations , '• J. Comp. Phys., Vol.57, 1985, pp. 327-360. [13] B. van Leer, "Towards the Ultimate Conservative Difference Scheme V : A Second-Order Sequel to Godunov's Methods," J. Comp. Phys., Vol. 32, 1979, pp. 101-13 6. [14] P. Colella and P.R. Woodward, "The Piecewise Parabolic Method (PPM) for Gasdynamics Simulations," J. Comp. Phys., Vol. 54, 1984, pp. 174-201. [15] P.L. Roe, "Approximate Riemann Solver, Parameter Vectors and Difference Scheme," J. Comp. Phys., Vol. 43, 1981, pp. 357-372. [16] P.L. Roe, "Numerical Algorithms for the Linear Wave Equations," Royal Aircraft Establishment Technical Report TR81047, 1981. [17] P.L. Roe and M.J. Baines, "Algorithms for Advection and Shock Problems," Proc. 4th GAMM Conference on Numerical Methods in Fluid Mechanics, H. Viviand, ed. , 1982. [18] A. Harten and J.M. Hyman, "A Self -Adjusting Grid for the Computation of Weak Solutions of Hyperbolic Conservation Laws," J. Comp. Phys., Vol. 50, 1983, pp. 235-269. [19] P.K. Sweby and M.J. Baines, "Convergence of Roe's Scheme for the General Nonlinear Scalar Wave Equations," J. Comp. Phys., Vol. 56, 1984, pp. 135-148. [20] S. Osher, "Riemann Solver, The Entropy Condition and Difference Approximations," SIAM J. Numer. Analy. Vol. 21, 1984, pp. 217-235. [21] S. Osher and S. Chakravarthy , "Very High Order Accurate TVD Schemes," The IMA Volumes in Mathematics and Its Applications, Vol. 2, Springer-Verlag, pp. 229-274, 1986. [22] H.C. Yee, "On the Implementation of a Class of Upwind Schemes for System of Hyperbolic Conservation Laws," NASA TM86839, Sept. 1985. [23] S.F. Davis, "TVD Finite Difference Schemes and Artificial Viscosity, " ICASE Report No. 84-20, 1984. [24] P.L. Roe, "Generalized Formulation of TVD Lax-Wendroff Schemes," ICASE Report NO. 84-53, Oct. 1984.

PAGE 147

136 [25] H.C. Yee, "Construction of Explicit and Implicit Symmetric TVD Schemes and Their Applications," NASA TM-86775 July 1985. Also J. Comp. Phys., Vol. 68, 1987, pp. 151-179. [26] C.J. Nietubicz, T.H. Pulliam and J.L, Steger, "Numerical Solutiion of the Azimuthal -Invariant Thin-Layer Navier-Stokes Equations," AIAA Paper 79-0010, Jan. 1979. [27] L.D. Kayser and F. Whiton, "Surface Pressure Measurement on a Boattailed Projectile Shape at Transonic Speed," U.S. Army Ballistic Research Laboratory, Aberdeen Proving Ground, Maryland, ARBRL-MR-03161, 1982. [28] P.D. Lax, Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves . Regional Conference Series in Applied Mathematics, 11, SIAM, Philadephia, 1973. [29] O.A. Oleinik, "Discontinuous Solutions of Non-linear Differential Equations," Uspekhi Math. Nauk. , Vol. 12, No. 3, 1957, pp. 3-73. ( Amer. Math. Soc. Transl., Ser. 2, Vol. 26, pp. 95-172.) [30] J.B. Goodman and R.J. LeVeque, "On the Accuracy of Stable Schemes for 2-D Scalar Conservation Laws," New York University Report, May 1983. [31] P.L. Roe, "Some Contributions to the Modelling of Discontinuous Flows," Proceedings of the AMS-SIAM Summer Seminar on Large-scale Computation in Fluid Mechanics, June 27-July 8, 1983. Lectures in Applied Mathematics, 22, 1985. [32] B. van Leer, "Towards the Ultimate Conservative Difference Scheme, II. Monotonicity and Conservation Combined in a Second Order Scheme," J. Comp. Phys., Vol. 14, 1974. pp. 361-370. ' • [33] J.E. Fromm, "A Method for Reducing Dispersion in Convective Difference Schemes," J. Comp. Phys., Vol. 3, 1968. pp. 176-189. [34] P.K. Sweby, "High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws," SIAM J. Numer. Anal., Vol. 21, 1984, pp. 995-1011. [35] R.F. Warming and R.M. Beam, "Upwind Second Order Difference Schemes and Appliciations in Aerodynamics," AIAA J. Vol. 14, 1976, pp. 1241-1249. [36] T.H. Pulliam and J.L. Steger, "Implicit Finite Difference Simulations of Three-Dimensional Compressible Flows," AIAA J Vol. 18, 1980, pp. 159-167.

PAGE 148

137 [37] J.L. Steger, "Implicit Finite Difference Simulation of Flow about Arbitrary Geometries with Application to Airfoils," AIAA Paper 77-665, June 1977. [38] C.J. Nietubicz, "Navier-Stokes Computations for Conventional and Hollow Projectile Shapes at Transonic Velocities," AIAA Paper, 81-1262, 1981. [39] B.S. Baldwin and H. Lomax, "Thin Layer Approximation and Algebraic Model for Separated Turbulent Flows," AIAA Paper 78257, Jan. 1978. [40] J.L. Steger, C.J. Nietubicz and K.R. Heavey, "A General Curvilinear Grid Generation Program for Projectile Configurations," BRL Memorandum Report -MR-03 142, Oct. 1982. [41] J.L. Steger and D.S. Chaussee, "Generation of Body Fitted Coordinates Using Hyperbolic Differential Equations," SIAM J. Sci. Stat. Compt., Vol. 1, Dec. 1980, pp. 421-437 [42] H.C. Yee, "Upwind and Symmetric Shock-Capturing Schemes," NASA TM-89464, May 1987.

PAGE 149

BIOGRAPHICAL SKETCH Ming-Hsiung Chen was born on September 2, 1951, in Kaohsiung, Taiwain, Republic of China. After high school he went to Nation Chung-Hsing University, Taichung, Taiwan, where he received his bachelor's and Master's degrees in soil and water conservation in 1974 and 1976 respectively. Since then he spent two years in military service at Chungli, Taiwan. In the summer of 1981, he enrolled as a graduate student in the Department of Engineering Sciences at the University of Florida. In 1982 he was married to Ko-Hui Liu. 138

PAGE 150

I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of philosophy. a C.C. Hsu Chairman Professor of Aerospace Engineering, Mechanics & Engineering Science I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as dissertation for the degree of Doctor of Philosophy. L.E. Malvern Professor of Aerospace Engineering, Mechanics & Engineering Science I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. U.H. Kurzweg Professor of Aerospace Engineering, Mechanics & Engineering Science

PAGE 151

I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. ^Xz W. Shyy Associate Professor of Aerospace Engineering, Mechanics & Engineering Science I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. A.K. Varma Professor of Mathematics This dissertation was submitted to the Graduate Faculty of the College of Engineering and to the Graduate School and was accepted as partial fulfillment of the requirements for the degree of Doctor of Philosophy. December, 1988 // t J, A A IJJLJ ft. & Dean, College of Engineering Dean, Graduate School

PAGE 152

UNIVERSITY OF FLORIDA 3 1262 08554 2925