
Citation 
 Permanent Link:
 http://ufdc.ufl.edu/AA00004795/00001
Material Information
 Title:
 On total variation diminishing schemes for transonic turbulent flow computation
 Creator:
 Chen, MingHsiung, 1951
 Publication Date:
 1988
 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 NavierStokes equations ( lcsh ) Turbulence ( lcsh )
 Genre:
 bibliography ( marcgt )
nonfiction ( marcgt )
Notes
 Thesis:
 Thesis (Ph. D.)University of Florida, 1988.
 Bibliography:
 Includes bibliographical references.
 General Note:
 Typescript.
 General Note:
 Vita.
 Statement of Responsibility:
 by MingHsiung 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 nonprofit 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 UTF8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID EOKPRX157_XJMVF0 INGEST_TIME 20111121T19: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 MINGHSIUNG 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. ChenChi 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 2D 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 UrbanaChampaign. Finally, the deepest gratitude is expressed to his wife KoHui 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 THREEDIMENSIONAL SYMMETRIC AND UPWIND TVD SCHEMES IN GENERALIZED COORDINATES 63 5.1 The 3D Equation in Generalized Coordinates. .. 63 5.2 Algorithm of 3D Schemes in Generalized Coordinates 72 5.3 A Conservative Linearized ADI Form for SteadyState Appliciation 75 VI . THIN LAYER NAVIERSTOKES SOLVERS 78 6 . 1 ThinLayer 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 SpatialVarying 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: NAVIERSTOKES EQUATIONS AND THIN LAYER APPROXIMATION 124 APPENDIX B: BEAMWARMING 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 "Ctype" grid for the SOCBT projectile with sting 93 The surface pressure coefficients of BeamWarming 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) BeamWarming 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 BeamWarming 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 BeamWarming 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 MINGHSIUNG CHEN December, 1988 Chairman : ChenChi Hsu Major Department : Aerospace Engineering, Mechanics & Engineering Science Recently a number of new techniques for constructing nonlinear, secondorder accurate, highresolution 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 twodimensional planar flows to
PAGE 11
threedimensional flows. On the application side, a thinlayer NavierStokes flow code based on the conventional BeamWarming 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 BeamWarming scheme. Furthermore, although the computing cost of the TVD schemes is higher than the BeamWarming 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 highspeed 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 wellknown 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 highfreguency 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 fourthorder explicit dissipation is used to control nonlinear instabilities, whereas the secondorder implicit dissipation is builtin to relax the stability bound from the explicit fourthorder dissipation. All these numerical dissipation terms contain adjustable parameters. For inviscid steadystate 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 BeamWarming scheme and other schemes .
PAGE 14
3 In order to overcome these difficulties, a number of new techniques for construction of nonlinear, secondorder accurate, highresolution 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 gasdynamic equations in characteristic form was proposed by Courant et al. [7]. Their procedure, sometimes called the CIR (CourantIssacsonRees) 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 firstorder accurate. Godunov's scheme In 1959, Godunov [8] proposed the idea of advancing the solution to the next timelevel 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 semiinfinite states (u = u L for x < 0, u= u R for x > ) For onedimensional 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 monotonicitypreserving 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 twostep 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 firstorder accurate TVD scheme. He then added a limited antidiffusive flux to a first order accurate entropysatisfying scheme to obtain the secondorder accurate TVD scheme. This technique is referred as the 'modified flux approach'. The modified flux is chosen so that the scheme is secondorder at a smooth region and will switch itself into firstorder 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 oneparameter family of implicit secondorder TVD schemes. More recently Harten's TVD scheme was modified and generalized by Yee [12] and implemented to solve the twodimensional 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 secondorder accuracy in Godunov's scheme by replacing the piecewiseconstant initial data of the Riemann Problem with piecewiselinear initial data. The slope of the piecewiselinear 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 shockcapturing 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 LaxWendroff 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 oneparameter family of secondorder 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 steadystate solutions are independent of the timestep. 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 threepoint central difference code ( e.g. BeamWarming 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 timedependent axisymmetric NavierStokes 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 timestep procedure [4] is also incorporated into this new version of the code to speedup the convergence rate for steadystate flow computations. The transonic turbulent flow over the axisymmetric secantogivecylinderboattailed (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 BeamWarming 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 twodimensions to threedimensions with generalized curvilinear coordinates. In Chapter VI, a thinlayer NavierStokes 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 variationdiminishing (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 RankineHugoniot 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, uc, Huc) 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 ik+l/ u ik+2/# u i+k) h il/2 = h ( u ik/ ui_ k+1 ,..., Ui+kx) 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 finitedifference 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 firstorder 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 oneparameter family of fivepoint 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(l0) (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 (l0)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) (Ru)i = ui A(l0) (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 (10) (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 i2/ u iif 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(l0)C i+1/2 >
PAGE 35
24 + .+ (2.27b) C i+1/2 + C i+1/2 = A (10) (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 secondorder shockcapturing schemes.
PAGE 36
CHAPTER III UPWIND AND SYMMETRIC TVD SCHEMES In this chapter, the description of nonlinear, upwind and symmetric, explicit and implicit, secondorder 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 firstorder 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 secondorder accurate TVD schemes by adding limited antidiffusive flux to a firstorder 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 FirstOrde r Accurate TVD Schemp Consider the onedimensional 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 threepoint 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 CourantFriedrichsLewy (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 nonphysical (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 CFLcondition is limited to use of a small time step. Therefore, in order to obtain steadystate solutions, it usually requires long computing time. This penalty is particulary severe for high Reynolds number NavierStokes 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 oneparameter 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 (10) (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 firstorder accurate in space. When 6 = 1/2, the scheme is secondorder 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 finitedifference operators, defined as (3.13a) (L.u)i = ui + A0(h i+1/2 hi_ 1/2 ) (3.13b) (Ru)i = ui A(l0)(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 CFLlike 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 SecondOrder Accurate TVD Scheme In order to achieve a secondorder accurate scheme which satisfies the TVD conditions, Harten modified the numerical flux fi of a firstorder 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 firstorder error. Therefore, Harten 's method is sometimes referred to as the modifiedflux approach. The modified flux is chosen so that the scheme is secondorder
PAGE 42
31 in smooth regions and firstorder at points of extrema. After introducing the modified numerical flux, the fivepoint, secondorder 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+l9i)/(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 steadystate calculation. One can rewrite g^ as gi = minmod(ai +1 /2Ai +1 / 2 u, <7ii/2Ml/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 socalled 'flux limiter' function for the control of unwanted oscillations in the numerical scheme. After careful examination of Harten's modifiedflux 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 firstorder 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 socalled 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 = <(7l)[ 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 twodimensional 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 xdirection, 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 ydirection. With the above notation, a oneparameter 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 RoeSweby's secondorder TVD scheme can be interpreted as a LaxWendroff 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 nonupwind 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 RoeDavis 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 secondorder 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(aii/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 nonconservative 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 steadystate application, one way to simplify Eq.(3.30) is to use a spatially firstorder 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 twodimensional implicit operators, such as Eg. (3.29a), is still very expensive. One can simplify the solution process by introducing an approximate factorization of the twodimensional operator into two onedimensional 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 Eii/ 2 ,j] " AY tFi,j+i/ 2 " Fi,jl/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 ydirection. Each sweep requires the solution of a linear system involving a block tridiagonal matrix which can be solved very efficiently by block LUD (lowerupper decomposition) solver.
PAGE 57
CHAPTER IV FLUX LIMITERS As already mentioned, firstorder accurate schemes, in general, suffer from numerical diffusion, white the classical secondorder accurate schemes such as the LaxWendroff 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 antidiffusive 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 threepoint 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 secondorder LaxWendroff scheme could be constructed by adding an antidiffusive flux to the firstorder upwind scheme, that is n+1 n n n (4.9 ) ui = Ui ^Ai_ 1/2 u a_[ H v(lv)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 LaxWendroff scheme produces spurious oscillations near the discontinuities of the solution, while the firstorder upwind scheme does not. By limiting the amount of the antidiffusive flux in Eq. (4.9), it is possible to obtain the secondorder accurate TVD scheme. The concept of this antidiffusive 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 twostep procedure. The antidiffusive 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 (lv)Ai +1 / 2 U ] The ^i is chosen to be a function of the ratio of consecutive cell gradients of A i1/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 + %(lv)[*(ri)/ri 4>{rii)\)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 LaxWendroff scheme, and if 4> (r) = r, Eq. (4.12) reduces to the secondorder 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 secondorder accurate is imposed by requiring that the amount of limited flux added to the first upwind scheme be a weighted averaged of the secondorder centered LaxWendroff scheme and the secondorder 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 secondorder schemes to produce a secondorder TVD scheme. He combined a monotonicity preserving but nonconservative version of the LaxWendroff scheme and the secondorder 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 WarmingBeam 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_ i1/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 secondorder TVD region of Figure 4.3. In [17], Roe proposed a secondorder 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 firstorder 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 nonphysical solution. Therefore using different limiters for different characteristic fields could be a good compromise. For the onedimensional 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 multidimensional system of conservation laws. In 1984, Davis [23] interpreted the upwind TVD finite difference schemes analyzed by Sweby [34] as LaxWendroff 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(lv)Ai +1 / 2 u %*j (11* I) (iQii/2)Aii/ 2 u h\u\ (1M) (lQi+l/2)Ai+l/2U. Here the first line represents the original LaxWendrof 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 iV2 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^(li/)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)Qil/2/*il/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 THREEDIMENSIONAL SYMMETRIC AND UPWIND TVD SCHEMES IN GENERALIZED COORDINATES An original derivation of the implicit, secondorder 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 3D Euler Equation in Generalized Coordinates The inviscid threedimensional 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 = (71) [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)Wf 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 threedimensional 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, Uk^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, Vk^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, Wk 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\lk z p k z U+kyp u+k x c uk x c k*v+k z p /CyV k z Vk x p V+kyC VkyC k x wkyp kyW+k x p ic z w W+k z C wk z c k x q 2 /2+ P (k Z VicyW) k y q 2 /2+ p (k x wk z u) k z q 2 /2+ P (ic y uk x v) H+C0 Ec'e (5.9) 1 R* = K x( 1 ~ lD l)~ Â«X b 2 u K z /p + ~Ky/p+ Â«x b 2 (Â« Z VKyW)/p " X b 2 v Â« x b 2 w /Cy(lbl)~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 /(7D + q 2 /2, bl = *>2 (q 2 /2), t>2 = (7l)/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 ith 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 ith 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 ith 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(yl) + 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 3D TVD schemes in generalized Coordinates With the above notation, a oneparameter 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 ^,kl/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 secondorder 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 steadystate 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 Â«il/2 + l<*i+l/2 Â»il/2l)/ i i (Â«i+l/2 + a il/2) / (5.19c) gi = S. max[0, min(2 lÂ«i +1/2 I , S'ali/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 secondorder 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 il/2 + "i+1/2)] 5.3 A Conservative Linearized API Form for SteadyState Application A conservative linearized ADI form of Eq. (5.14) for steadystate 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 ji/ 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, jl/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 fourtorder 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 NAVIERSTOKES 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 Reynoldsaveraged NavierStokes equations in generalized curvlinear coordinates are simplified by using the thinlayer 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, noniterative approximate factorization finite difference scheme developed by Beam and Warming [3]. The axisymmtric thinlayer NavierStokes 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 NavierStokes 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 Azimuthalinvariant 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 = ji pvU pwU+e 2 p (e+p)UÂ£ t p pW pUW+f x p m 1 u f +m 2 f x 6 = J1 pvW pwW+f 2 p ; 3 = = ji m 1 w f +m 2 r z (e+p )Wf t pJ m 1 m 3 +m 2 (r x u+f A,
PAGE 92
81 (6.2c) ft = Jl and Â° pV[RÂ£(UÂ£ t )+R r (Wft)] pVR(Vi, 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 (7l)" 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 twolayer 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 backwarddifference 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 firstorder 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 fourthorder artificial dissipation term (6.4) D e = e e At Ji^A^V^) 2 + (A^Vj) 2 ] J, is added to the righthand side of Eq. (6.3) and the implicit secondorder 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 secondorder 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 (lowerupper 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 BeamWarming scheme and the TVD scheme is that the TVD scheme requires more sophisticated dissipation terms. The dissipation terms of the Beamwarming 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 Rii/ 2 *ii/2 + R k+l/2*k+l/2 " R kl/2*kl/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,ki/ 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 SecantOgiveCylinderBoattail (SOCBT) projectile shape is shown in Figure 6.1. Specifically, the projectile under study is 3caliber for secantogive 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 xz 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 secondorder extrapolation. Numerically this is done by using a threepoint 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 twopoint 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 noslip 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 SpatialVarying Time Step For steady state flow problems, it is desirable to achieve the converged solution as quickly as possibble. The convergence to a steadystate 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 BeamWarming 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 secantogivecylinderboattailed (SOCBT) projectile model is adopted as the test problem. An axisymmetric thinlayer NavierStokes 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 nonphysical 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 thinlayer NavierStokes 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 " Ctype" grid for a flow domain of 4 projectile lengths with special clustering on both the ogivecylinder and cylinderboattail 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 "Ctype" 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 BeamWarming scheme 90x60 "Ctype" grid I.M9MJ S.^mH LSXtIH (.MM lUÂ« Figure 7.2a The surface pressure coefficients of BeamWarming 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 thinlayer code, the parameters of usercontrolled 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 CRAYXMP/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 BeamWarming 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 BeamWarming 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 BeamWarming scheme, the proper choices of the free parameters such as q and e e are not always clear; a trialanderror 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 BeamWarming 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.wttn tr 2: cc o z CM ancB 0. JOItB MflCH NO. = 0.96 I I ' ' I I I L2MiO% O.fMfrfCM ITERATIONS V Figure 7.3 Convergence history of different schemes (a) BeamWarming 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 BeamWarming 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 BW 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) BW 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 : BeamWarming 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 BeamWarming 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 BeamWarming 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 BeamWarming scheme and symmetric TVD scheme. Overall, the BeamWarming 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 BeamWarming 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 Lonorm 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.mmm x/o Figure 7 . 5a The surface pressure coefficients of BeamWarming scheme for projectile with MÂ« = 0.96, and Re = 7.7xl0 5 using a 90x40 C grid. The L2norm residual of 2.59xl0" 6 can be reached in around 3000 time steps. Q MflCH NO. = 0.36 o ex peri netrial 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 steadystate 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.KKB .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 L2Norm 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.Mtcn o CM J 6.MCM 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. 10KB ' Â— ' 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, = 12, 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 rr 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.MCW O.MCM 0.WKB : 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 threedimensional generalized curvilinear coordinate systems has been presented. The derived schemes are then simplified and implemented into a timedependent axisymmetric thinlayer NavierStokes 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 BeamWarming 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 BeamWarming 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 NAVIERSTOKES EQUATIONS AND THIN LAYER APPROXIMATION The strong conservation law form of the threedimensional NavierStokes equations in Cartesian coordinates can be written in the flux vector form as (Al) dq 3E 3F 3G dt ax 3y dz 3 Ey a Fy + + ax ay aG v az where (A2) 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 (A3) 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 (7l)" 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 ( = (71) [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 NavierStokes Equations To enhance numerical accuracy and efficiency and to handle boundary conditions more easily, the NavierStokes equations can be transformed from Cartesian coordinates to general curvilinear coordinates where (A4) 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 (A5) 8r q + dÂ£ d 8 3 . (E A E V ) + (F F v ) + (6 6 V ) dr, 0, ar where q = Jl (A6) pU pW e E = Jl P U pUU + Â£ x p pvU + Â£ y p pWU + Â£ Z P (e+p)U Â£ t P
PAGE 138
I = Jl 127 pV puV + ry x p pW + r/yp pWV + Â»? z p (e+p)V r/ t P 6 = ji pW pUW + f x p pVW + f y p pWW + f z p (e+p)W rtP and (A7) 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 (A8) Â£ 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 heatflux vector in nondimensional form are given in Eq. (A3) . 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 Xy r ), (A9) 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 Â£ Jl = ThinLaver Approximation In highReynoldsnumber flows, the viscous effects are confined to a thin layer near noslip 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. (A5) simplifies to (A10) a T q +a^E +a^F +a r d = Re _1 a r Â§, where (All) 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 (A13) 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 BEAMWARMING IMPLICIT BLOCK ADI ALGORITHM AND ARTIFICIAL DISSIPATION The finitedifference algorithm due to Beam and Warming applied to Eq. (A10) results in the following approximate factorization : (Bl) 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 Reifijfinj D e , and where 6 is the centraldifference operator and A and v are forward and backwarddifference operators, and h = At corresponds to the firstorder timeaccurate Euler implicit scheme. DjJ^, Di  n and Di  jare the secondorder implicit and D e is the fourthorder 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 (B2) 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 + jlftn^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 = (B3) Â«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 KzCyep" 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 (B4) 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 (pi) 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/NavierStokes 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. (Bl) uses a constant coefficient dissipation model. In this approach, an explicit fourthorder dissipation (B5) D e q n = Â£e At ^[(V^) 2 + (V^A,,) 2 + (V r A r ) 2 ] Jqn, is added to the righthand side of the equation and an implicit secondorder dissipation (B6) Di k = ei jl 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. 217237. [2] R.W. MacCormack, "The Effect of viscosity in Hypervelocity Impact Cratering," AIAA Paper 69354, Cincinnati, Ohio, 1969. [3] R. Beam and R.F. Warming, "An Implicit Factored Scheme for the Compressible NavierStokes Equations," AIAA Paper 77645, 1977. also AIAA J. Vol. 16, 1978, pp. 393401. [4] T.H. Pulliam and J. Steger, "Recent Improvements in Efficiency, Accuracy and Convergence for Implicit Approximate Factorization Algorithms," AIAA Paper 850360, 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. 357393. [6] A. Harten, J.M. Hyman and P.D. Lax, "On FiniteDifference Approximations and Entropy Conditions for Shocks," Comm. Pure Appl. Math., Vol. 29, 1976, pp. 297322. [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. 271306. Also Cornell Aeronautical Lab. translation. [9] J. P. Boris and D.L. Books, "FluxCorrect Transport I SHASTA : A Fluid Transport Algorithm That Works," J. Comp. Phys., Vol. 11, 1973, pp. 3869. [10] S.T. Zalesak, "Fully Multidimensional Flux Correct Transport Algorithms for Fluid," J. Comp. Phys., Vol. 31, 1979, pp. 335362. [11] A. Harten, "On a Class of High Resolution TotalVariationStable Finite Difference Schemes," SIAM J. Numer. Anal. Vol. 21, 1984, pp. 123. 134
PAGE 146
135 [12] H.C. Yee, R.F. Warming and A. Harten, "Implicit Total Variation Diminishing (TVD) Schemes for SteadyState Calculations , 'Â• J. Comp. Phys., Vol.57, 1985, pp. 327360. [13] B. van Leer, "Towards the Ultimate Conservative Difference Scheme V : A SecondOrder Sequel to Godunov's Methods," J. Comp. Phys., Vol. 32, 1979, pp. 10113 6. [14] P. Colella and P.R. Woodward, "The Piecewise Parabolic Method (PPM) for Gasdynamics Simulations," J. Comp. Phys., Vol. 54, 1984, pp. 174201. [15] P.L. Roe, "Approximate Riemann Solver, Parameter Vectors and Difference Scheme," J. Comp. Phys., Vol. 43, 1981, pp. 357372. [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. 235269. [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. 135148. [20] S. Osher, "Riemann Solver, The Entropy Condition and Difference Approximations," SIAM J. Numer. Analy. Vol. 21, 1984, pp. 217235. [21] S. Osher and S. Chakravarthy , "Very High Order Accurate TVD Schemes," The IMA Volumes in Mathematics and Its Applications, Vol. 2, SpringerVerlag, pp. 229274, 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. 8420, 1984. [24] P.L. Roe, "Generalized Formulation of TVD LaxWendroff Schemes," ICASE Report NO. 8453, Oct. 1984.
PAGE 147
136 [25] H.C. Yee, "Construction of Explicit and Implicit Symmetric TVD Schemes and Their Applications," NASA TM86775 July 1985. Also J. Comp. Phys., Vol. 68, 1987, pp. 151179. [26] C.J. Nietubicz, T.H. Pulliam and J.L, Steger, "Numerical Solutiion of the Azimuthal Invariant ThinLayer NavierStokes Equations," AIAA Paper 790010, 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, ARBRLMR03161, 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 Nonlinear Differential Equations," Uspekhi Math. Nauk. , Vol. 12, No. 3, 1957, pp. 373. ( Amer. Math. Soc. Transl., Ser. 2, Vol. 26, pp. 95172.) [30] J.B. Goodman and R.J. LeVeque, "On the Accuracy of Stable Schemes for 2D Scalar Conservation Laws," New York University Report, May 1983. [31] P.L. Roe, "Some Contributions to the Modelling of Discontinuous Flows," Proceedings of the AMSSIAM Summer Seminar on Largescale Computation in Fluid Mechanics, June 27July 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. 361370. ' Â• [33] J.E. Fromm, "A Method for Reducing Dispersion in Convective Difference Schemes," J. Comp. Phys., Vol. 3, 1968. pp. 176189. [34] P.K. Sweby, "High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws," SIAM J. Numer. Anal., Vol. 21, 1984, pp. 9951011. [35] R.F. Warming and R.M. Beam, "Upwind Second Order Difference Schemes and Appliciations in Aerodynamics," AIAA J. Vol. 14, 1976, pp. 12411249. [36] T.H. Pulliam and J.L. Steger, "Implicit Finite Difference Simulations of ThreeDimensional Compressible Flows," AIAA J Vol. 18, 1980, pp. 159167.
PAGE 148
137 [37] J.L. Steger, "Implicit Finite Difference Simulation of Flow about Arbitrary Geometries with Application to Airfoils," AIAA Paper 77665, June 1977. [38] C.J. Nietubicz, "NavierStokes Computations for Conventional and Hollow Projectile Shapes at Transonic Velocities," AIAA Paper, 811262, 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 MR03 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. 421437 [42] H.C. Yee, "Upwind and Symmetric ShockCapturing Schemes," NASA TM89464, May 1987.
PAGE 149
BIOGRAPHICAL SKETCH MingHsiung Chen was born on September 2, 1951, in Kaohsiung, Taiwain, Republic of China. After high school he went to Nation ChungHsing 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 KoHui 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

