Citation |

- Permanent Link:
- http://ufdc.ufl.edu/AA00025703/00001
## Material Information- Title:
- A self-adaptive computational method for transonic turbulent projectile aerodynamics
- Creator:
- Reed, Christopher William, 1960-
- Publication Date:
- 1987
- Language:
- English
- Physical Description:
- vii, 108 leaves : ill. ; 28 cm.
## Subjects- Subjects / Keywords:
- Approximation ( jstor )
Coordinate systems ( jstor ) Geometric lines ( jstor ) Grid generation ( jstor ) Heating ( jstor ) Mathematical variables ( jstor ) Orthogonality ( jstor ) Projectiles ( jstor ) Transonic flow ( jstor ) TVD schemes ( jstor ) Aerodynamics, Transonic ( lcsh ) Curvilinear coordinates ( lcsh ) Dissertations, Academic -- Engineering Sciences -- UF Engineering Sciences thesis Ph. D Fluid dynamics ( lcsh ) Numerical grid generation (Numerical analysis) ( lcsh ) Projectiles -- Aerodynamics ( lcsh ) - Genre:
- bibliography ( marcgt )
theses ( marcgt ) non-fiction ( marcgt )
## Notes- Thesis:
- Thesis (Ph. D.)--University of Florida, 1987.
- Bibliography:
- Includes bibliographical references (leaves 105-107).
- Additional Physical Form:
- Also available online.
- General Note:
- Typescript.
- General Note:
- Vita.
- Statement of Responsibility:
- by Christopher William Reed.
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- Copyright [name of dissertation author]. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
- Resource Identifier:
- 021625499 ( ALEPH )
18145334 ( OCLC )
## UFDC Membership |

Downloads |

## This item has the following downloads: |

Full Text |

A SELF-ADAPrIVE CURyTIONAL METHOD FOR TRANSONIC TURBULENT PROJECTILE AERODYNAMICS BY CHRISTOPHER WILLIAM REED 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 1987 To Shelley and Drew ACKNOWLEDGEMENTS The support of Dr. Chen-Chi Hsu has been most encouraging during the course of this work. His guidance and willingness to provide it are appreciated. In giving me the freedom to experiment, discover, fail and succeed I have learned as much about research as I have about the topic. For this I am grateful. My appreciation is extended as well to the committee members, Dr. Leadon, Dr. Kurzweg, Dr. Mikolaitis and Dr. Wilson. This work was partially supported by the AFOSR. Computer time was provided by NASA Ames, NSF and the Pittsburgh Supercomputer Center. TABLE OF CONTENTS ACKNOWLEDGEMENTS . . . . . . . . . . . ABSTRACT . . . . . . . . . . . . . CHAPTER I INTRODUCTION . . . . . . . . . . II THE ADAPTIVE GRID GENERATION EQUATIONS . . . . 11.1 The Curvilinear Coordinate System . . . . 11.2 A Variational Approach . . . . . . 11.3 Inversion of the Grid Generation Equations . III NUMERICAL SOLUTION OF THE ADAPTIVE GRID GENERATION EQUATIONS . . . . . . . . 111.1 The Newton-Raphson Iterative Method . . . 111.2 Modified Solution Procedure ............ III.3 The Adaptive Grid Generation Code ......... IV CONTROL FUNCTIONS . . . . . . . . . IV. 1 The Basic Form . . . . . . . . IV.2 Sxmoothing . . . . . . . . . . IV. 3 Other Control Functions . . . . . . V FLOW SOLVERS . . . . . . . . . . V.1 The Governing Equations ................. V.2 Solution Algorithms for the Governing Equations V.3 Examples on Fixed Grid Networks . . . . VI PRELIMINARY TESTS AND MODIFICATIONS . . . . . VI.1 The General Procedure . . . . . . . VI. 2 Orthogonality . . . . . . . . . VI. 3 Curvature . . . . . . . . . VI.4 Internal Grid Boundaries . . . . . . VII THE SELF-ADAPTIVE COMPIRTATIONAL METHOD . . . . iii vi 1 7 7 12 17 22 22 26 * 37 39 39 41 42 44 44 46 51 61 61 * 62 62 69 76 D J VIII RESULTS AND DISCUSSION . . . . . . . . . 81 VIII. 1 Choices for the Control Functions . . . . 81 VIII.2 Results Using the Self-Adaptive Computational Technique: The Beam and Warming Scheme . . .. .84 VIII.3 Results Using the Self-Adaptive Ccaputational Technique: The TVD Scheme . . . . . .. .90 VIII.4 Applications to the Projectile Base Flow Problem 97 IX CONCIJDING REMARKS . . . . . . . . . .. 102 APPENDIX . . . . . . . . . . . . .. 104 REFERENCES . . . . . . . . . . . . .. 105 BIOGRAHIICAL SKETCH . . . . . . . . . . .. 108 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy A SELF-ADAPTIVE CtXPTATIONAL METHOD FOR TRANSONIC TURBJIENT PROJECTILE AERODYNAMICS By Christopher William Reed August, 1987 Chairman: Chen-Chi Hsu Major Department: Engineering Sciences An adaptive grid generation procedure is developed for viscous flow problems. The equations governing the adaptation are based on a variational statement resulting in a set of elliptic equations in which adaptation can occur independently in each coordinate direction. The method allows for explicit control of adaptation and orthogonality while grid smoothness is implicit in the elliptic equations. The adaptive grid generation equations retain a simple relationship between the control functions and the grid point spacing, are capable of efficiently providing the extremely refined mesh in the boundary layer regions and can maintain a smooth grid network in complex geometries. A self-adaptive computational technique has been developed by coupling the adaptive grid generation procedure with a thin layer Navier-Stokes code and a TVD code. In solving an axisymmetric turbulent transonic projectile flow problem, the grid network was adapted to the pressure distribution in the streamwise direction and to the velocity distribution in the direction normal to the projectile surface. Results obtained for Mach 0.91, 0.96 and 1.2 indicate that the procedure can provide good adapted grid networks provided good choices are made for the control functions. vii CHAPTER I TDRJCTION The study of boundary fitted grid systems has become an important part of ccmpitational fluid dynamics and, in fact, their development has been cited as a major pacing item in the progress of computational fluid dynamics [1]. The use of such systems makes the application of finite difference techniques versatile and effective in solving complex fluid dynamic problems in arbitrary domains. Boundary fitted grid generation can be viewed as the development of a general curvilinear coordinate system in which the coordinates coincide with the boundaries on the edges of the physical domain. This feature of the boundary fitted grid systems facilitates the application of finite difference approximations at the boundaries and eliminates the need for interpolation between grid points near the boundary. This is particularly advantageous when the boundaries are highly curved or have slope discontinuities since the use of interpolation near these boundaries may have a significant effect on the solution [2]. In the domain interior, the intersection of each family of coordinate lines denotes each grid point and therefore defines the computational mesh. The distribution of the coordinate lines in the physical domain may 2 be concentrated to provide high resolution in certain areas but will always form an equally spaced rectilinear coordinate system in the transformed plane. Thus, once the governing equations are transformed onto the curvilinear coordinates the finite difference algorithms used to solve these equations may be developed on the equally spaced rectilinear coordinate system which is inherent in the transformed plane. Certain characteristics of the general curvilinear coordinate system have been shown to affect the accuracy of the finite difference approximations in the transformed plane [3]. These include orthogonality of intersecting coordinate lines and the smoothness of the grid spacing. It is also known that a refined mesh is needed in regions of large physical gradients to obtain accurate approximations. In fact, poor resolution in these high gradient regions can be detrimental to solution accuracy and to the convergence of finite difference algorithms for fluid flow problems [4]. Thus, the development of a good adapted grid network is important for the successful application of finite difference methods to solve complex fluid dynamic problems. The self-adaptive grid generation technique has become an important area in computational fluid dynamics since it has been shown to provide good grid networks for the complex flow fields occurring in transonic and supersonic flows [5,6,7]. The use of boundary fitted curvilinear coordinate systems with transformed governing equations leads naturally to the concept of solution adaptive grids. As constraints of computer storage and CPU time prohibit uniform mesh refinement throughout the entire domain, the coordinate spacing in the physical domain is varied to increase resolution in only those regions in which the solution is changing rapidly. For simple flow problems, when the position of important solution gradients is known, good adapted grids can be obtained with conventional techniques. However, in more complex problems the position of these important regions is not known a priori and then the development of a proper adaptive grid network is difficult. Solution adaptive grid generation addresses this problem by continuously updating the grid network as the solution evolves such that the important physical gradients are sufficiently resolved as they develop. It is thus an attempt to efficiently reduce the truncation error by only refining the mesh in those regions where the error may be large. Also, since the grid characteristics of smoothness and orthogonality also affect the truncation error, a good adaptive grid generation technique should include the enhancement of these characteristics as well. The general approach of most adaptive grid generation schemes is based on minimization techniques. A measure of each desired grid characteristic is defined, and the governing equations for the grid are obtained by minimizing the integral of these measures over the domain. In many instances adaptation is required in only one coordinate direction. Dwyer et al. [8], for instance, have developed an equidistribution law to control the grid spacing along one coordinate direction. The spacing was adapted to a flame front in a two-dimensional combustion problem and the improved resolution eliminated spacial oscillations in the flame front. Gnoffo [9] has developed another one dimensional adaptive scheme in which tension springs are assumed to connect the grid points along a coordinate. The spring constants became large when increased resolution was required, and the points clustered in those regions in an attempt to minimize the potential energy of the spring system. The application of this technique to supersonic flow over a Viking Aeroshell successfully 4 resolved a separated shear layer by adapting to the velocity gradient. One dimensional adaptation can be extended to higher dimensions by successive adaptation in each coordinate direction. Nakahashi and Deiwert [7] have used the spring analogy in such an approach and added torsional springs that connect intersecting coordinate lines to control orthogonality. They solved a transonic viscous airfoil problem in which the grid was adapted to the density gradient in the streamwise direction to resolve shocks and adapted to the velocity gradient normal to the airfoil surface to resolve the boundary layer. The one-dimensional approach to multi-dimension adaptation has the advantage of efficiency since the grid can usually be obtained using a marching solution algorithm. Another advantage is the independence of adaptation in each coordinate direction. However, it is difficult to maintain smoothness in such approaches, and highly skewed grids may result [10]. Rai and Anderson [11] have used a different approach for two dimensional adaptation. In their scheme each grid point induces a velocity, either repulsive or attractive, on every other point. By choosing some flow gradient to indicate the need for point clustering, each point with a gradient larger than the average gradient attracts points, and those with values below the average repel points. By summing the induced velocities from each point, the velocity of each point is determined and can be integrated over a time step to determine the new point location. In applying this scheme to the two dimensional linearized viscous Burger's equation, they adapted the grid to the derivative of the dependent variable in each direction and showed that such an approach will increase solution accuracy. 5 Saltzman and Brackbill [5] have developed an adaptive grid generation scheme based on a variational approach. A functional is defined to measure each grid characteristic and the minimization of these functionals results in a set of partial differential equations that govern the grid point spacing. With the proper choice of parameters the equations are elliptic, which guarantees a one-to-one mapping and results in a smooth grid. They solved an inviscid supersonic flow problem by adapting the grid cell size to a function of the pressure gradient to capture shocks. However, it is not possible in this approach to adapt the grid spacing independently in each coordinate direction, a feature that is prohibitive if the spacing in different coordinate directions vary by a large amount. The adaptive grid generation technique presented here is similar in approach to that of Brackbill and Saltzman but is developed for applications to viscous transonic flow problems. The solution of these problems usually contains shock structures aligned with one coordinate direction and boundary shear layers parallel to a streamwise coordinate. Also, the grid point spacing can vary by orders of magnitude along different coordinates and it is thus necessary to adapt the grid independently in each coordinate direction. Therefore, a functional is defined for independent adaptation in each coordinate direction such that the resulting equations are elliptic. The use of elliptic equations in grid generation offer some important advantages that motivate their use. They constitute a boundary value problem and can thus be used in any domain. They guarantee a one to one mapping and also, grid smKothness is inherent in the equations. In addition to the equations governing the grid, there are two other topics that require attention when developing 6 an adaptive grid generation technique. First is the definition of the control functions, which are the part of the grid generation equations that dictate the need for adaptation. They are dependent upon the solution and thus provide the link between the grid generation and the governing equations. Also, a method is needed to account for the effect of the moving grid on the solution. Since the solution is defined at each point, any movement will alter the solution. The proposed adaptive grid generation technique is used to adapt the grid for a transonic axisynmmetric projectile flow problem which is solved using both an implicit factorized algorithm and a TVD scheme for the thin layer Navier Stokes equations. The constraint of axisynmmetric flow reduces the domain to two independent spatial variables, thus the grid generation equations are developed in two dimensions. The technique, however, can be readily extended to three dimensions. CHAPrER II THE ADAPTIVE GRID GENERATION EQUATIONS 11.1 The Curvilinear Coordinate System The use of boundary fitted curvilinear coordinate systems in solving equations on arbitrarily shaped domains was originally motivated due to problems in applying finite difference approximations at curved boundaries. Because curved boundaries do not in general coincide with the rectilinear Cartesian coordinates, the grid points defined along these coordinates do not usually lie on the bounding curve. It therefore becomes difficult to accurately represent the boundaries [2] and special procedures must be used in those regions [12]. A solution to this problem, which has become an effective and versatile approach, is to map the arbitrarily shaped physical domain (x,y) into a rectangular computational domain ( Y7). Such a mapping for a typical projectile flow domain is shown in Figure 2.1. As indicated, the curved C and 17 coordinate lines coincident with the boundaries in the physical domain map onto the edges of the rectangular computational domain. The boundary grid points, being defined by the coordinate intersections, are coincident with the boundary in the computational plane and thus no special treatment for the boundary conditions is required. Another advantage is that virtually any set of nonuniformally spaced curvilinear coordinates in the physical domain map into an equally spaced regular z--fm (1.JM) j (,.,1M; Ij.1 ________ _________________________ _________i ___ ___ ____ ________ __ __ j__ [___' t.1) F Figure 2.1 Curvilinear coordinate system '1M.! ) 4 9 mesh in the computational domain. Furthermore, once a successful finite difference algorithm is developed on the computational domain, any physical domain may be mapped into this region and solved with the same numerical algorithm. The task of numerical grid generation is to define the transformation between the Cartesian and computational coordinate systems which can be viewed as the development of a curvilinear coordinate system in the physical domain. The general mapping between the two coordinate systems is C = (x,y) (2.1.1) 1= 7 (x,y) Owing to the discrete nature of finite difference approximations, it is sufficient to define a finite set of coordinate lines, the intersections of which define the grid points in the physical domain. Thus for each grid point (x,y) in the physical domain, there are two curvilinear coordinates ( t) that intersect at that point and designate its location. It is convenient to view this mapping in the computational domain, where for each pair (C,t7) there corresponds a point (x,y) in the physical domain. This correspondence can be conveniently expressed as the inverse mapping x = x(C,v) (2.1.2) y = y(,?I) 10 It is advantageous to designate the curvilinear coordinates with integer values. Finite difference expressions in the computational plane are simplified since A and At become unity and the mapping of eqs. (2.1.2) will appear as a two dimensional array with integer subscripts, creating a naturally ordered way to store the grid points. Using the inverse mapping, the location of any coordinate line in the physical domain can be determined simply as the set of points obtained by holding one of the arguments in the mapping constant. The integer indicies i and j are used here to designated each grid point. Thus the ith point along a coordinate is i and the jth point along an n coordinate is 1j. For simplicity the notation (xi,j,yij) is used to represent the point (x(i,qj),Y(i,?j)) and the total number of points in each direction and r are IM and JM, respectively. For the projectile problem considered, the and i coordinate lines will always run in the streamwise and normal directions as indicated in Figure 2.1. Once a mapping is defined, in a discrete sense, the governing equations can be obtained in the computational plane by transforming them onto the curvilinear coordinate system. Derivatives with respect to the Cartesian coordinates can be expressed using eqs. (2.1.1) and the chain rule as [ax a { } I { } (2.1.3) | a |y a |7 Iay Jy [Qr where x, y, nx, and ?7y are the metric coefficients of the transformation. They appear as constants in the transformed equations and depend solely on the mapping of eqs. (2.1.1). Substitution of eqs. 11 (2.1.3) into the governing equations makes and n the independent spatial variables and now a numerical algorithm can be developed using finite difference approximations on the equally spaced computational domain. Since the metric coefficients depend only on the mapping, virtually any domain may be mapped into the computational domain and solved using the numerical algorithm. The metric coefficients can be evaluated for a numerically generated transformation by considering the metrics of the inverse transformation. The relationship between the metric and inverse metrics is, for two dimensions, {x y i = 1 Y -Y J (2.1.4) J y 17y -x X where J is the Jacobian of the transformation, J = x y7 xy Y (2.1.5) The derivation of these equations can be found in reference [13]. Since and r are the independent variables for the inverse metrics, the inverse metrics can be evaluated numerically in the computational domain and subsequently used to determine the metric coefficients. The primary disadvantage of using boundary fitted curvilinear coordinate systems is the appearance of additional truncation error terms in the finite difference approximations on the computational plane. Mastin [3] has shown, for instance, that the truncation error T for a second order finite difference approximation of the first derivative is, in the ccxrputational domain T = -x2 f xxx fxx (2.1.6) T=6x fxx 2 The first term is the usual truncation error term containing the grid point spacing x and the third derivative of the function. The second term is due to the varied spacing along the $ coordinate lines in the physical domain. It is thus possible to increase the truncation error if the coordinates spacing is not smooth. They have also shown that truncation error for a first derivative varies inversely as the sine of the angle between the curvilinear coordinates and therefore a lack of orthogonality can also increase the truncation error. Orthogonality is also important near solid boundaries for turbulent flow problems because turbulence models often are based on information in the normal direction. In the context of a coordinate transformation, the adaptive gridding can be viewed as an attempt to reduce the truncation error by adjusting the transformation metrics. This is accomplished by judiciously defining the curvilinear coordinates such that they vary smoothly, are nearly orthogonal, and concentrate in regions where the solution is changing rapidly. 11.2 A Variational Approach The choice of a method to generate the curvilinear coordinates is arbitrary in that it has no dependence on the partial differential equations that are to be solved. Of the many methods available, however, the variational approach is attractive for purposes of adaptive grid generation because it can provide a simple, explicit means of incorporating desired properties into a grid generation scheme and can 13 result in a set of elliptic equations which govern the grid. To develop an adaptive grid generation scheme using variational techniques, a measure of each grid characteristic is defined and its integral is taken over the physical domain. If the measures are defined such that they measure a deviation from the desired property then the coordinates and ,7 that minimize the integral will yield a grid with the desired characteristics. As a minimization problem, the Euler-Lagrange equations can be applied to obtain a set of partial differential equations, the solution of which defines the curvilinear coordinates. The first of three functionals used here is Ip = J v8dxdy (2.2.1) which measures the adaptation of the grid spacing in the coordinate direction to the control function P. When P is small, the quantity V must also be small in order to minimize the integral; a small value of v7 corresponds to large grid spacing. Consequently, when P is large the grid spacing along the coordinate will be small. Similarly, the integral Iq= J Q Vndxdy (2.2.2) is a measure of the adaptation of the grid spacing in the P direction to the control function Q. The third functional Io = (v.vi1)2 dxdy (2.2.3) 14 measures the orthogonality of the and q coordinates. As will be shown later, grid smoothness can be controlled through the definitions of the control functions. These three integrals can be combined into one total functional IT IT { V'V + Vn.Vn + (V7.V?7)2 } dxdy (2.2.4) I P Q where x is a parameter that weighs the relative importance of adaptation and orthogonality. Before applying the Euler-Lagrange equations for and q, it is beneficial to scale each term in the equations so that they remain of the same order of magnitude. An ordering analysis shows that the first term in eq. (2.2.4) is of order (C2/P), the second is of order (C2/Q) and the orthogonality term is of order (C4/L2) where C is a computational length scale and L is a physical length scale. It is ccmnnon practice to use highly stretched grids in transonic flow problems to provide extremely small grid spacing in the viscous sublayer region. The spacing along the normal coordinate can increase by five orders of magnitude, consequently the terms in eq. (2.2.4) may vary in a similar way. It is therefore beneficial to scale the terms locally rather than globally since a global scale cannot reflect the size of the term throughout the domain. The local scales Ap, Aq, and Ao used here are Ap=PL Aq=QL Ao=JL (2.2.5) where PL and QL are the local values of the control functions and JL is the local Jacobian and is of order (L2/C2). Although they vary throughout the domain they are considered constants during the application of the 15 Euler-Lagrange equations. This violates the variational principal, but it has been shown [13] that the use of local scaling will produce better grids when there are large changes in the grid spacing. After substituting in the scales, the total functional IT becomes IT= J ( PL v + QLVn .V + AJL (V'.Vq)2 )dxdy P Q (2.2.6) After applying the Euler-Lagrange equations for and P to the total functional IT while holding the scales constant, the following set of partial differential equations is obtained: xx + yy + Vp.V + A(n2 xx+2qxqy~xy+?yy + (2Cxx+Cyqy) Y7xx + (t x y+ y x) xy+(x x+2yry) yy)=O (2.2.7) 'lxx + 'lyy + V.Vn + ((2xtx+yty)xx+( x~y+qy~x)xy + Q 2 (Qx'x+2Cy'y) xy + x'2xx+2x~ytjxy + C217yy)=0 The subscript L has been dropped from the local scales in the differential equations since they are no longer needed to indicate their assumed invariance. The and ri coordinates that satisfy these equations will yield a grid with the desired properties defined in the functionals of eqs. (2.2.1), (2.2.2) and (2.2.3). It should be noted that there was no explicit functional to maintain grid smoothness. There are two types of smoothness that should be differentiated. One type, which can be controlled by the control functions, is smoothness of the grid point 16 spacing along a coordinate line. If the control functions change rapidly along a coordinate, the grid spacing will consequently change rapidly. This smoothing thus can be governed through the definition of the control functions. Another type of smoothness, which is inherent in the elliptic equations, is that between the grid point spacing of adjacent coordinates. If the grid point distribution along two adjacent coordinates varies, the elliptic equations will tend to dampen the variation. Equations (2.2.7) govern the grid coordinates and n and will be elliptic as long as the parameter A is not too large. They constitute a boundary value problem and it is therefore necessary to define the coordinates along the boundaries. It is also necessary to adapt the boundary points in a consistent manner with the interior points so that the boundary points remain aligned with the interior points. For this purpose, a one dimensional functional analogous to equation (2.2.1) or (2.2.2) is defined to measure adaptation along a boundary coordinate, IS= f s ds (2.2.8) where s is the arclength along the boundary coordinate. There is no consideration for orthogonality and scaling is unnecessary. Applying the one dimensional Euler-Lagrange equation yields a second order differential equation, Sss-Cs Ps/P=O (2.2.9) Note that a similar equation can be obtained for boundary adaptation along a r coordinate by replacing in eq. (2.2.9) with ? and P with Q. This one dimensional equation, together with eqs. (2.2.7), are the 17 equations that govern the grid coordinates in the physical domain. 11.3 Inversion of the Grid Generation Equations In order to solve the grid generation equations, they are transformed onto the curvilinear coordinate system. This is done by inverting eqs. (2.2.7) and (2.2.9) to make x, y, and s the dependent variables. The curvilinear coordinates then become the independent variables and the grid generation equations are solved on the computational domain which is equivalent to defining the inverse mapping of eqs. (2.1.2). All the advantages of using the computational domain to solve the governing equations discussed previously will be realized here as well. There are two approaches to inverting the equations. One is to transform the functional IT onto the computational domain and then apply the Euler-Lagrange equations for x and y. Another approach, used here is to derive the corresponding expressions on the computational domain for each term in eqs. (2.2.7) and (2.2.9) and then to substitute them into the equations. The first term in eq. (2.2.7) can be written using eq. (2.1.3) as aa xx = ( x s-- + "X -- ) X (2.3.1) Using eq. (2.1.4) to transform the metrics this expression becomes Y,18 yC a Y'i Cxx =( - ) (2.3.2) J a9 J -9i J and after differentiating, the relationship is found to be -Y__ 2 x 2 2 1XX- J 2(y2x-2ygyx -y+2x)+-(y~y -2ygyEy7+yy2 (2.3.3) Similar expressions can be derived for each of the other second derivative terms in eq. (2.2.7), they are listed in Appendix A. All the first derivative terms, of course, can be transformed using eqs. (2.1.4). By substituting these expressions into eqs. (2.2.7) and collecting like terms the equations for adaptive grid generation become, in the cciputational domain, alxE + a2x n + a3x, + blyE + b2Yn + b3Yq = S1 (2.3.4) clxC + C2X, + c3Y, + dly + d2Yn + d3Ynq = S2 in which a, = -y7(c( + A'#2) 2A'yeaf a2 = 2y (# + A'f-y) + A'y (4,2+J2) a3 = -yi (I + A'72) 2A'y af3 (2.3.5) bI = x,7 (a + A',2) 21yct( b2 = -2x,, ( + A'7-y) A'x (4#2+J2) b3 = x,(7-y + A'-y2) + 2A'y73fi CI = yE(a + A'a2) + 2A'ynap C2 = -2y (6 + A'), ) A'y, (4,32+J2) C3 = y$(-y + A',32) + 2A'yn7p di = -xa (Q + A'a2) + 2A'yaf3 d2 = 2xE(62 + A'ap) + A'y(462+J2) d3 = -Y (7 + A'\2) + 21'yr7 19 S= x2 + Yb2 Y = x2 + y2 = x Y + xy A' = A/J and PCJa P7J S1= +-- P p (2.3.6) Q Jf QJ7 S2 -= + Q Q It is important to note that in inverting the equations the derivative of P in the Y direction appears even though P was intended to control the spacing in the direction. This occurs because the contravariant base vector v actually indicates the spacing in the direction normal to lines of constant rather than in the direction. If the grid is orthogonal, the two directions, v1 and V, are equivalent and the ) derivative will not have any influence since 6 will vanish. It has been found in numerical experiments that eliminating this term from the equations has a negligible effect and consequently it is dropped from the equations. The source terms S1 and S2 in eqs. (2.3.5) then become S1 =-J PJa s=^ (2.3.7) QtJ-y S2 -=Q Q 20 The one dimensional equation for the boundary adaptation, eq. (2.2.9), is also transformed in a similar manner. The second derivative term can be shown, using the chain rule, to transform as Sss - 3-- SC3 (2.3.8) Substituting eq. (2.3.7) into eq. (2.2.9) yields the transformed one dimensional equation for adaptation along C boundaries (2.3.9) s6 + s PC/P =0 and for adaptation along boundaries coincident with 7 coordinates snn + SnQn/Q = 0 (2.3.10) These equations can be solved directly for the grid point spacing along a coordinate. For instance, integrating eq (2.3.10) once yields the grid point spacing S = C/Q (2.3.11) where c is a constant of integration and can be shown to be ST C --- f da? JQ (2.3.12) where ST is the total arclength along the coordinate in the physical domain. However, since the solution in the interior will require an iterative algorithm, the differential form, eqs. (2.3.9) and (2.3.10) are 21 retained and solved iteratively so that the boundary points remain consistent with the interior points as the grid evolves. CHAPTER III NUMERICAL SOI=ION OF THE ADAPTIVE GRID GENERATION EQUATIONS III.l The Newton-Raphson Iterative Method Equations (2.3.3), (2.3.9) and (2.3.10) constitute a set of nonlinear elliptic partial differential equations. In order to solve these equations numerically they are approximated with second order central difference expressions which results in a set of coupled algebraic equations for each grid point. The finite difference approximation for the equations at grid point, (xij, Yi,j) become al(xi+l,j b, (yi+l, j cl (xi+l,j dl (yi+l, j -2xi,j + xi-l,j) + a3(xi,j+1 -2xi,j + xi,j-I) + -2yij +Yi-l,j) + b3(Yi,j+l -2Yi,j + 2Yi,j-l)-FI=O0 (3.1.1la) -2xi,j + Xi-i,j) + c3(xij+l -2xij +xi,j-l) + -2Yi,j +Yi-l,j) + d3(Yi,j+l -2Yi,j +Yi,j-l)-F2=0 (3.1.1b) PJcx F1 = 2 aT1 b2T2 (3.1.Pc) (3.1i.IC) Q J7 Q?2 1-- a i d2T2 F2 a2T1 d2T2 23 T1 = (Xi+l,j+l + Xi-l,j-l1 Xi+l,j-l xi-l,j+l)/4 (3.1.id) T2 = (Yi+l,j+l + Yi-l,j-l Yi+l,j-l Yi-lj+l)/4 where the source terms and the coefficients all contain first derivatives which when approximated with central differences do not contain the point xij, Yi,j explicitly. This set of equations is solved using a Newton- Raphson iterative scheme. The initial guess for the grid point locations will not satisfy the finite difference representation of the equations and therefore a non-zero residue will exist (Rl)i,j = eq (3.1.la) (3.1.2) 2 (R2)i,j = eq (3.1.1b) where 2 indicates the iteration number. An estimate of the residue at the next iteration can be obtained by expanding the residue in a Taylor series about the current residue and retaining only the linear terms, 2+1 2 a(Rl) i,j a(RI)i,j (Rl)i,j = (Rl)i,j + Axi,j -- + Ayi,j -- (3.1.3) S2 2 2+1 1 a(R2)i,j 8(R2)i,j (R2)i,j = (R2)i,j + AXi,j3 + Ayij, Where Ax and y are the changes in the grid point locations where AX and Ay are the changes in the grid point location 24 AYi,j =Yi,j-Yi,j (3.1.4) The derivatives of the residues with respect to the variables Xji,j and Yi,j can be obtained by differentiating eqs. (3.1.1) and are (R )i ,j axilj a (R2)i ,j aXilj An approximate residue at the = -2(al+a3), = -2(c1+c3), value of Axij and Ayi, j next iteration to vanish I(Rz) ] (R2i,)j .(R2) l,j a (Rl)-l,j axi,j a (R2)1,j axi,j a (Ri,j a(R2)i,j ayi,j = -2(bl+b3) (3.1.5) = -2(dl+d3) can be obtained if we require the a(RI)1,j 8Yi,j a (R2)i ,j lYi,j {Axij) Ayi,j (3.1.6) Substituting eqs. (3.1.5) into eqs. (3.1.6) Axij and Ayij yields f Axi,j 1 2 r (dl+d3) I AYi,j J D -(c+c3) and solving for -(bl+b3) 1*(Rl)i,,j (al+a3) (R2) il,j (3.1.7) where D = 4((al+a3)(dl+d3)-(cl+C3)(dl+d3)) (3.1.8) 25 Knowing Axij and Ayij, the new point location can be calculated by 2+1 2+i solving eqs. (3.1.4) for xij and Yij. This same procedure can be applied to the one dimensional equations for adaptation, eqs. (2.3.9) and (2.3.10). The residue for eq. (2.3.9) is (Rs)i = si+1 2si + Si-1 + sP/P and the equation for the new point location becomes s2+l-2+^=2 1 1i/2i (3.1.9) (3.1.10) (3.1.11) The numerical algorithm consists of updating the value of each point on the boundary with eq. (3.1.11) and of sweeping through the domain and updating each interior point using eq. (3.1.7). The criterion used here to indicate convergence is Rmax < 6 (3.1.12) where Rmax is the maximum residue in the interior weighted by the Jacobian AX2 j + AYiJ )1/2 Rmax= max( '--' )1/2 iJ (3.1.13) Currently, the criterion 6=0.005 is used. 26 In the original iterative procedure, the most recently updated values for each grid point location were used in updating the remaining points and thus the procedure resembles a Gauss-Seidel iteration. This approach, however, is not suitable for vector processing and thus cannot benefit from the increased processing speed. Therefore, another approach is used in which all the grid points along a coordinate line are updated in a Jacobi iterative fashion. The new grid point values along that coordinate are then used in updating the points along the next coordinate line. Thus the technique can be viewed as a half Jacobi, half Gauss-Seidel iteration. This technique, which is suitable for vector processing, reduced the computer processing time for each iteration by 40 percent. Furthermore, only a few more iterations were required for the same convergence criterion in comparison to the Gauss-Seidel iteration, thus the effective speed up in processing time is approximately 40 percent. III.2 Modified Solution Procedure As pointed out earlier, the grids used for viscous transonic flow problems contain highly refined grid spacing in the direction normal to the surface (P? direction) in order to resolve the viscous sublayer. This spacing results in grid cells with aspect ratios, up to 105 in the present applications, near solid boundaries. Iterative solution procedures for elliptic equations become inefficient for such large aspect ratios since the motion in one direction will be limited by the small spacing in the other direction. The effect of the aspect ratio on an iterative process can be ascertained by considering the simple grid cell shown in Figure 3.1 which for convenience is rectangular and aligned with the Cartesian (Lb) LX wu)__ ) (1,0 Figure 3.1 High aspect ratio cell f *1~1 ________ 1 _____ -- -. ______ - __ __ ii- irni - - ___ -~ I ~ i iii--, Figure 3.2 Initial grid network y,1 4, O, b ) (2,b) (2.0) 1,0,01 28 coordinates. For this example the aspect ratio AR is equal to (1/b). Assume for purposes of demonstration that the two terms (Q /Q) and (P/P) are non-zero at the center point and will therefore cause the point to move. We can obtain the changes in the point location by evaluating eqs. (3.1.7) for this simple case. Many terms for this example are zero xY=0 xCn yC =0 y =0 y 17 =0 x =0 ye,=o (3.2.1) yn=0 yr =0 The coefficients become al=b3 a2=0 a3=b bl=0 b2=0 b3=0 cl=0 c2=0 c3=0 dl=b2 d2=0 d3=l (3.2.2) and upon substitution into eqs. (3.1.7) yield P b Ax'j 2P 1 l+b Q1 1 Ay, ( ) -2Q l+b (3.2.3) Writing these in terms of the aspect ratio AR, its effect on the grid point motion becomes clear. If we compare the solution for a square cell , i.e. AR=1, 29 AX (.) )xi 2P 2 (3.2.4) Aij=- ( ) Yi,' 2Q 2 with that for a cell of large aspect ratio AR >> 1, Pe( 1 Axjij 2P l+AR (3.2.5) oi7 AR Ay, =1 ( AR ) Yij 2Q 1+AR it becomes evident that the reduction in the grid point motion in the x direction for large AR is proportional to (1/AR). For aspect ratios of the order (105) this result becomes prohibitive, rendering the iterative procedure extremely inefficient. An example has been formulated to demonstrate this deficiency. In an initial grid network, as shown in Figure 3.2, the points along the coordinates are clustered in one region and the points along each v coordinate are clustered near the bottom surface to form the highly stretched spacing typical of grid networks used in transonic flow calculations. The spacing at the surface is approximately 10-5 times that in the upper region of the grid. In the adaptive grid generation equations, the control function P is set equal to a constant so that the points will redistribute themselves with equal spacing along the coordinates. Using a technique described in Chapter IV, the control function Q is defined so that the highly stretched spacing along each 30 coordinate will be retained in the solution to the adaptive grid generation equations. Using the Newton-Raphson iterative solution procedure, the solution is advanced 20 iterations which results in the grid shown in Figure 3.3. As can be seen, the point spacing along the coordinates is approximately equal for the coordinates sufficiently above the surface. However, in the region of the high aspect ratio cells, the movement of the points in the direction is retarded by the extremely small spacing in the i direction. A modification has been made to eq. (3.1.11) which updates the points along the bounding coordinate in formulating this example. The movement of points using the one dimensional equation is independent of the interior points and thus is not restricted by the high aspect ratio cells. This creates an inconsistency at the boundary since the points adjacent to the boundary are constrained by the high aspect ratio cells. Equation (3.1.11) therefore has been modified as A SA 1 ASi = (Rsi)i (-- ) (3.2.6) 1+AR. where AR* is the aspect ratio of the cell adjacent to the ith point along the boundary coordinate. This change only effects the rate of convergence of the iterative solution procedure for the boundary points so that their movement will remain consistent with that of the interior points adjacent to the boundary during the convergence process. In order to alleviate the deficiency of the above solution procedure, a technique has been developed which is based on solving for a reduced grid in which many of the points in the 7 coordinate direction have been ",I \ \ \ I I Ill/ \L/ \ II I! Figure 3.3 Grid network after 20 iterations -- __ -_ - 7M - 7= Figure 3.4 Reduced aspect ratio cell I 32 removed to decrease the maximum cell aspect ratio. Figure 3.4 shows a diagram of a typical grid for viscous flow problems containing large aspect ratio cells near a solid boundary. A reduced grid can be formed by removing all the coordinate lines denoted with the dashed lines, resulting in a grid with fewer points along the r7 coordinate lines which are designated by the intersection of the solid lines. Currently, enough points are removed along the q coordinate lines such that the minimum spacing between coordinate lines is greater than 0.04. The typical spacing along the boundary is of order (0.1) which results in cell aspect ratios of order (1). Upon convergence of the solution for the reduced grid, all the points removed are reinserted along r coordinate lines using the same one dimensional equation that controls the spacing along the boundaries. It is necessary in this process to adjust the control function Q at the remaining points to account for the spacing requirements of the points removed. A simple procedure for this purpose can be derived by using the equation for one dimensional adaptation, eq. (2.3.10). The procedure will also be adequate for the two dimensional adaptive grid generation equations. If we write explicitly the finite difference approximation for eq (2.3.10), sj+-j+1j-l Qj+I-Qj-l Sj+l-2sj + sj-i + (- 2 )( 2 )/Qj=0 (3.2.7) we can regroup the terms in the form Qj+I-Qj-l QJ + (--2- -) (Sj+l-Sj)=(sj-sj-l) (- Qj+---- ) (3.2.8) Q (Qj+I-Qj-l 2 33 which, after using forward difference expressions to indicate the spacing between two adjacent points can be written ASj = ASj-1 Cj Qj+i-Qj -i QJ + ( 2 ) Cj= 2) Qj+I-Qj-I Qj ( --- -)Q 2 (3.2.9) (3.2.10) Now, as shown in Figure 3.5, if the ith point is removed it will be necessary to adjust the control function such that the new relationship is satisfied ASj+i = AS*-1 C3-I (3.2.11) where Cj-1 is the modified control function and AS3-i is the increased grid point spacing * ASjl = ASj-|- + ASj (3.2.12) The correct form of Cj-l can be obtained in the following derivation ASj+1 = ASj_- Cj-i ASj Cj+1 = (ASj-l+ASj) Cj-1 Asj-1 Cj Cj+1 = ASj-1 (1+Cj) Cj-1 (3.2.13) cj511 c- + 1+ _1 I+Cj J- 34 Once C3-_ is calculated using eq. (3.2.13), the modified value of Qj- can be obtained by solving eq. (3.2.10) for Qj-i l+C*-I Q*- = (- ) (Qj+l-Qj-2) (3.2.14) l-C*-1 This procedure is applied for each point along an j coordinate that is removed. It should be noted that this procedure does not affect the final grid, the points remaining in the reduced grid will obtain the same location as if all the points were used in the iterative solution procedure. The removed points are, once the solution for the reduced grid is obtained, inserted back along the r7 coordinates lines. This is done by approximating the r? coordinate with straight line segments between the points of the reduced grid and then reinserting the removed points along each segment with spacing dictated by eqs. (2.3.10) using the original control function Q. If necessary, the use of a higher order approximation to the r coordinate lines, such as a cubic spline, can be used to obtain a better representation of the n coordinate lines. In using this procedure with the two dimensional grid generation equations, the solution may vary somewhat since the equations are coupled. Also, the effects of orthogonality may alter the results in some areas. In order to show the advantage of this procedure, the initial grid network of the previous example is advanced 20 iterations using the reduced grid technique. Figure 3.6 shows the initial reduced grid in which the points responsible for the high aspect ratio cells are removed. After modifying the control function Q, the solution to the adaptive grid generation equations is advanced 20 iterations resulting in the grid shown in Figure 3.7. As can be seen, the points along all coordinates Figure 3.5 Corresponding grid point spacing ___Figure 3.6 Initial reduced grid network_ __ Figure 3.6 Initial reduced grid network I Figure 3.7 Reduced grid network after 20 iterations Figure 3.7 Reduced grid network after 20 iterations Figure 3.8 Final grid network after reinserting grid points 37 are approximately equally spaced and the spacing along the r coordinate lines in the initial reduced grid is retained which indicates that the modified control function is correct. Upon inserting the removed points along the ri coordinates using eq. (2.3.12), the final grid network, Figure 3.8, is obtained, which is the same grid that would have been obtained in Figure 3.3 after a large number of iterations. III.3 The Adaptive Grid Generation Code A computer code base on the iterative solution procedure described in the preceding section has been developed and is referred to here as the Adaptive Grid code. As the solution process is iterative, an initial guess for the grid point locations is required and the flow variables for which the grid network is to be adapted are needed for calculation of the control functions. These are considered as input to the Adaptive Grid code. The initial guess for the grid network can be a grid network obtained by any conventional grid generation technique or a previously adapted grid network. Values for the flow variables at the grid points are obtained using the flow solvers which are described in Chapter V. The sequence of calculations in the Adaptive Grid code proceeds as follows: 1. Input an initial grid network and the flow variables at each grid point. 2. Calculate the control functions P and Q. 3. Remove points to form the reduced grid network and calculate the modified control function Q*. 4. Solve the grid generation equations for the reduced grid network using P and Q* . 5. Reinsert the removed points along the n coordinate lines using the control function Q which then results in the adapted grid network. 38 The methods used to calculate the control functions from the solution variables is described in Chapter IV. CHAPTER IV CONTROL FUNCTIONS IV. 1 The Basic Form It is the role of the control functions in the equations governing grid adaptation to dictate the need for increased resolution and, as indicated in section 11.1, they should be judicially chosen so that the resulting grid resolution reduces what would otherwise be large truncation error. For instance, Pierson and Kutler [14] chose for the control function, when solving a one dimensional inviscid Burger's equation, the third derivative of the dependent variable which is also the leading truncation error term in their finite difference approximation to Burger's equation. Their one dimensional adaptive grid scheme was thus an attempt to minimize a measure of the truncation error over the domain. In more complex problems, however, with more dependent variables and in higher dimensions, the truncation error expressions are complex and cumbersome to calculate. Therefore, in practice, an understanding of the physical problem and experience guide the choices for the control functions. The general control function considered here is P = 1 + Jpf (4.1.1) Q = 1 + 7qg 40 where 7p and -Yq are parameters and f and g are sane derivatives of the dependent variables scaled to range between zero and one. In the following developments, the terms W,-y and h will be used to indicate either P,7p and f or Q,-q and g. By evaluating numerically the solution to the one dimensional equation for adaptation, the following expression equating the control function to the grid point spacing between points (i) and (i+l) along a coordinate line is obtained, 1 ST1 ASi = +-yhi (4.1.2) 1 j l+-yhj After evaluating this expression for the maximum and minimum spacing corresponding to h = 1 and h = 0 respectively, and forming their ratio (ASi)max ( m 1 + -y (4.1.3) (ASi)mim it can be seen that the specification of 7 dictates the ratio of maximum and minimum spacing along a coordinate line. The parameter 7-y also influences the smoothness of the resulting grid point distribution. If it is chosen as zero, equal spacing as indicated by eq. (4.1.2) would result; as 7 is increased, the spacing must change more rapidly to obtain the varied spacing. Following the work of Nakahashi and Deiwert [7], another parameter a has been introduced to obtain more control over the spacing. The general control function now becomes W = 1 + 7ha (4.1.4) 41 By prescribing the minimum spacing and 7-y, eq (4.1.2) can be written as 1 ST (415) (ASi)min (l+jh ) -1+ (4.1.5) in which a is the unknown quantity and can be determined by a root finding technique. Currently, the Newton-Raphson finding method is employed. Thus, both the minimum and maximum grid point spacing along a coordinate can be specified by specifying the two quantities, (Asi)min and -Y. It should be noted that the prescribed values may not be realized in practice due to the coupling effects of the two dimensional adaptive grid generation equations. Also, the constraint of orthogonality may alter the results in same areas. IV. 2 Smoothing The choices for the functions f and g will of course depend on the problem being solved. However, the functions will in general, contain derivatives of the solution which are approximated numerically and it is therefore useful to smooth the control functions to eliminate any spurious data. Also, since smoothing is most effective when the function changes most rapidly, it will help reduce any rapid variations in the control function and subsequently smooth the grid point spacing. The method of smoothing used here, which is similar to that used in reference [5] is Wi+ij+Wi-ij ^J i=2,IM-I Wij= WiJ + 2 Wij j=I,JM (4.2.1) f Wi,j+l+Wi-lj I i=lIM Wij= WilJ + -W j=2,JM-I 'j=W~ 2 1,I where the value of W used here is 0.4. The control function at a point is replaced by a weighted average of itself and the two adjacent points along a coordinate line. Currently, two consecutive sweeps of this algorithm are made over the domain for each of the control functions, P and Q. In this algorithm there is no accounting for the varying spacing between points and thus this smoothing is in the computational plane. IV. 3 Other Control Functions In many instances a certain point distribution in one coordinate direction is known to yield good results and it would be advantageous to maintain this prescribed distribution in that direction while allowing adaptation to the solution in the other directions. For instance, in some of the viscous transonic flow problems solved here, it is known that a distribution in the direction normal to the surface based on an exponential function will yield good results provided enough points are used. A procedure to incorporate specified point distributions into the adaptive grid scheme is easily developed by using eq. (2.3.10). Consider some given grid point distribution s in the n coordinate direction which is to be reproduced by the adaptive grid generation equations 43 S = S(7) (4.3.1) After solving for the control function, eq. (2.3.10) becomes 1 Qi = (4.3.2) (s,)j The control function then can be determined from the prescribed distribution s by using central difference approximations to evaluate S,. As noted previously, the resulting distribution may vary somewhat due to the coupling of the two dimensional adaptive grid equations and effects of orthogonality. This approach, which is used in some of the cases solved here to specify the spacing in the normal coordinate direction, is similar to the technique of Thompson et al. [15] to obtain boundary point distribution in the domain interior. The equation for the grid point spacing used to obtain the arclength distribution is Asj = As*(l+e)n (4.3.3) where e is determined so that ST matches the total length of the coordinate line (see reference 16), As* is a specified spacing at the projectile surface, and n is equal to the index j. The actual choices for the functions h and g in eqs. (4.1.4) are discussed in Chapter VIII. CHAPTER V FLOW SOLVERS V.1 The Governing Equations The governing equations for viscous transonic flow are the Navier- Stokes equations and an equation of state. The Navier Stokes equations, which represent the conservation of mass, momentum and energy form a set of parabolic/hyperbolic partial differential equations when the time dependent terms are retained. To solve these equations numerically on an arbitrary domain the governing equations are scaled and then transformed onto the curvilinear coordinates. For problems in three dimensions, three curvilinear coordinates are required, denoted here as , q and . However, for bodies of revolution, such as the axisymmetric projectiles at zero angle of attack considered here, the solution is assumed invariant in the circumferential direction, designated as the direction, and the equations can be reduced since only two independent spatial variables are required. A further reduction in the equations ccmes from the thin layer approximation. In the thin layer approximation all the viscous derivatives in the direction along a solid surface (the direction in the present applications) are dropped from the equations. Due to restrictions on computer storage and CPU time, it is not usually feasible to resolve the viscous terms in both the directions normal to the surface and along the surface. For high Reynolds number flows, such 45 as the transonic flow considered here, large velocity gradients are expected to occur normal to the surface and most of the grid points are used to resolve these gradients. As larger spacing is then used along the surface and the viscous derivatives in that direction are expected to be small in comparison, they are dropped from the equations. The conservative form of the transformed thin layer governing equations for axisynmetric flow, written in the compact vector form are atq + a8 + a ? + t =Rela (5.1.1) where AP q p pu q=J* ( pv | PW c SpU A pUU+ xp , E=J pvU+ yp P pwU+ zp S(e+p) U ] A S, F=J 3 SpV p uV+n xP ( pvV+yP pwV+7I zp S(e+p)V 0 6=J* ( -pW2y p/y [ pW(yU+yV) I 0 (5.1.2) 0 uRu, + (u/3)T'7x uRv 7 + (u/3)T'y uRw, + (u/3)TVz A( +(u2+v +w2) +uPr1(-1)-i(a2)2 +(u/3) ( 7xu+, 2yv+t2xw)T and R = 'x2 + ny2 + nz2 T = ?,u + n yv, + nxW, A S = J 46 Here, p is the density, p is the pressure and E is the total specific energy. The velocities u,v and w correspond to the x, y and z directions, respectively, as indicated in Figure 5.1 and U, V, and W are the contravariant velocity components corresponding to the , t7 and directions respectively. The coefficient of viscosity is u, Pr is the Prandtl number, Re is the Reynolds number, -y is the ratio of specific heats and a is the local speed of sound. The details of the transformation and application of the thin layer approximation to obtain eqs. (5.1.1) are available in reference [13]. The assumption of a calorically and thermally perfect gas was made in closing the system of equations. The coefficient of viscosity pa is considered variable and it is also necessary to account for the effects of turbulence for the present applications. Under the assumption of invariance, only a two dimensional grid network containing and 7 coordinates in a plane through the axis of symmetry will be required. It should be noted, however, that the metric coefficients represent the three dimensional coordinate transformation and are therefore different from the expressions of eqs. (2.1.4) and the Jacobian J* now represents the grid cell volume. V.2 Solution AlQorithms for the Governing Equations The implicit spatially factored numerical scheme of Beam and Warming [17] is used to solve the transformed governing equations. The scheme is implemented for three dimensional problems in a code due to Steger and Pulliam [18]; the code used here is a modified version which was developed by Nietubicz et al. [19] for the efficient calculation of axisymmetric projectile problems and is referred to here as the Thin Figure 5.1 Coordinate systems for projectile configuration C Figure 5.2 Boundary configuration for projectile with sting 48 Layer code. The code employs Sutherland's law [20] to relate the viscosity and heat conductivity to the temperature and the effects of turbulence are modeled using the eddy viscosity model of Baldwin and Lomax [21]. The numerical algorithm is a tine-marching scheme in which steady state solutions are obtained as the time asymptotic solution of the unsteady equations. Spatial derivatives are approximated with central differences and a forward difference approximation is used for the unsteady terms. The time marching scheme of the finite difference approximation to the governing equations are written in delta form as (I + h 6. eID1) (I + h6n hRelJkJ1 ejD2A (5.2.1) = h(dC + 8, Re1 6nA 6) ED3 where A and B are the Jacobian matrices a A A a B (5.2.2) aq a and M is obtained as the time linearization of S. Here, 6 is the central difference operator. Details of this scheme can be found in Warming and Beam's [17] work, as well as in that of others. The terms D1, and D2 are implicit dissipation operators and D3 is an explicit dissipation term D1 = J-1 A v7 J D2 = J-1 An VI J (5.2.3) D3 = J-1 ((As V7)2 + A, V7)2) j 49 where D and v are forward and backward difference operators, respectively. Since the scheme is not naturally dissipative these terms are added to the finite difference scheme to prevent nonlinear instabilities and high frequency spatial oscillations. The constant coefficients fI and EE used here are 4At and 2At respectively which are consistent with values used in other applications of this code [19]. In general, these coefficients should be chosen large enough to dampen the instabilities but kept small enough to prevent any degradation of solution accuracy. In Chapter VII, same examples using the Thin Layer code with the Adaptive Grid code show that there are problems in obtaining a converged solution for the thin layer equations when the grid is adapted. Therefore another numerical scheme for solving the thin layer equations is also used. This scheme is the TVD scheme developed by Yee [22] for the multidimensional Navier-Stokes equations. The theory and development of TVD schemes is beyond the scope of current discussions. However, for purposes here it is sufficient to view the TVD scheme as an improvement in the artificial damping terms in the Beam and Warming scheme. In fact, the TVD scheme developed by Yee for multidimensions can be implemented into the existing Thin Layer code based on the Beam and Warming scheme by modifying only the added dissipation terms. The original Thin-Layer code has been modified to provide an option for the TVD scheme [23] and is referred to here as the TVD code. These 'sophisticated' dissipation terms switch the scheme from second order to first order accuracy at points of extrema and result in a more dissipative scheme in those regions. It should be pointed out that the primary purpose of the current research is to develop an adaptive grid generation technique. The TVD scheme is used 50 here because of problems occurring when using the Beam and Warming scheme in conjunction with the adaptive grid generation technique. These problems are discussed more fully in Chapter VII. The first flow problem considered is a the axisynmetric projectile flow problem with an attached sting to eliminate the base flow region. In both the Thin Layer code and the TVD code the boundary conditions are implemented explicitly by updating the flow variables at the boundary after each time step. Referring to Figure 5.2 which shows the basic C- type grid configuration for the axisymmetric projectile flow problem, there are four different sets of conditions applied at the four bounding coordinates. Along the outer boundary, line AB, free stream conditions are assumed. Thus the outer boundary is required to be sufficiently far from the projectile so that this assumption is valid. Along the upstream line of symmetry, line AD, the variables are obtained by requiring symmetry about the x axis. Numerically this is done by setting a second order forward difference approximation to the first derivative along the coordinate for each variable equal to zero and solving for the values at the points on the axis in terms of known values in the domain. Along the dowmstream boundary, line BC, the variables are obtained by extrapolation along the coordinate lines. However, if the free stream velocity is subsonic the pressure is fixed. This is done to prevent oscillations in the pressure distribution that would otherwise occur [24]. Along the projectile surface, line CD, the no slip condition holds for viscous flow, thus the velocity components are zero. The density is obtained by zero order extrapolation along the 7 coordinate lines and the pressure is obtained by satisfying the momentum equation along the surface. 51 Another problem considered is the projectile base flow problem. In Figure 5.3, which shows the basic 0-type grid configuration, the coordinate lines bend around the sharp corner at the base and end on a downstream axis of symmetry which is now an q coordinate. The existing codes are used for this problem by simply modifying some of the boundary conditions. The boundary conditions along the upstream axis of symmetry, line AD, and the projectile surface, line CD, are the same as those used in the projectile problem with sting. For the downstream axis of synnetry, line BC, the same conditions along the upstream axis are used; however, a backward finite difference formula is required. The freestream conditions are applied along the outer boundary along the line AA' and the downstream boundary conditions for the projectile with sting problem are applied along the portion of the outer boundary between A'B, except that the extrapolation is done along the r coordinate lines. V.3 Examples on Fixed Grid Networks In order to show the general characteristics of the Beam and Warming and TVD schemes, and for purposes of comparison with the solutions obtained with the adaptive grid generation technique, the two schemes were used to solve the projectile flow problem with sting on a fixed grid network. The fixed grid network was obtained using a method developed by Steger et al. [16] which is an effective technique for external flow problems. Two equations, one enforcing orthogonality and one specifying the grid cell volume, form a set of hyperbolic partial differential equations which are solved numerically to generate the grid network. In the marching solution algorithm, the initial grid points are specified along the projectile surface. The solution is marched outward in the 7 A A D C Figure 5.3 Boundary configuration for projectile base flow problem SOCBT 6.0 3.017 -- ----- 2.0 3.6-- ALL DIMENSIONS IN CALIBERS Figure 5.4 SOCBT projectile configuration 53 coordinate direction, defining each successive coordinate line such that the spacing between each coordinate line is consistent with the prescribed grid cell volume. At the bounding r coordinate lines, the upstream axis of symmetry and the downstream boundary of Figure 5.2, the coordinate lines are required to intersect the t coordinates at right angles. The projectile used in all the calculations is a 6-caliber, secant- ogive cylinder boattail configuration which is shown in Figure 5.4. There are experimentally determined pressure coefficients available at a series of points along the projectile surface for a range of Mach numbers in the transonic range [25]. The data will serve as a means of comparing the computed solutions. A grid network for this projectile problem with sting was generated using this technique with 70 points in the coordinate direction (IM=70) and 60 points in the 7 coordinate direction (JM=60). As shown in Figure 5.5, most of the points in the t? coordinate direction were clustered near the projectile surface to resolve the boundary layer. The spacing at the projectile surface in the i coordinate direction is 0.00002. The outer boundary is extended to four times the projectile length. In the streamwise direction, the points were clustered near the secant-ogive/cylinder and cylinder/boattail junctures in anticipation of the expansion waves that will occur in transonic flow. The first solution on this grid network is obtained with the Thin Layer code for Mach 0.96. The Reynolds number for this and all examples is 76,000. For all the results obtained, the time step used for a converged solution is 0.1. Figures 5.6 and 5.7 show the convergence process of the solution algorithm for the indicated time spans. It is quite obvious that the solution oscillates in the region containing the Figure 5.5 Initial fixed grid network for projectile with sting (70 x 60) MACH 0.96 0.2 0.0 -0.2 - -0.4- -0.6 - Figure 5.6 Cp 0.2- 0.0 - o experiment - computed Converging solution for the Beam and Warming scheme MACH 0.96 ;,, *..,,4 -0.2 / -0.4 / g i o experiment -0.6 1 7 computed Figure 5.7 Converged solution for the Beam and Warming scheme MACH 0.96 o experiment - computed Figure 5.8 R 1 10 10- Converged solution for the TVD scheme Beam and Warming 1 2 3 4 5 N (1000) Least squared residue 0.0 -0.4f -0.6 Figure 5.9 57 shocks. However, after approximately 5000 time steps the solution does converge and agrees well with the experimental data as indicated by the comparison of the experimental and computed pressure coefficient along the projectile surface. Figure 5.8 shows the results obtained using the TVD code. The solution also agrees well with the experimental data; moreover it converged within 1600 time steps. This difference in the convergence ofthe two schemes is shown in Figure 5.9 which plots the least square residue R, versus the number of time steps N. Here the slower, oscillatory convergence of the Beam and Warming scheme is evident while the TVD scheme appears to converge monotonically. The projectile base flow problem is also solved on a fixed grid using the TVD code. The grid network for this calculation was obtained using a modified version of the technique of Steger et al. [16] in which the boundary condition on the downstream bounding j coordinate line has been changed so that the coordinate lines intersect with the x axis at right angles. The grid network obtained with IM=70 and JM=60 is shown in Figure 5.10. As can be seen, the coordinate lines now bend around the sharp corner, and end at the downstream axis of symmetry forming an 0-type grid network. The computed pressure coefficient obtained for this problem using the TVD code is shown in Figure 5.11. Again, the solution converged in approximately 1600 time steps. It compares well, as before, over the secant-ogive and cylinder portions but then differs from the experimental results at the end of the boattail. No experimental data is available for the base region. There are a number of reasons that may cause the difference between the computed and experimental results near the corner. One reason may be the presence of a narrow sting used to support the experimental model that is not present in the computational 58 configuration. The difference could also be due to poor resolution in the projectile corner region or an inadequate turbulence model near the flow separation region. Figure 5.10 Initial fixed grid network for projectile base flow problem (70 x 60) MACH 0.96 (j Cjx2e'ment S Or0) U i ed MACH 0.96 - b a Ii 0.0 Figure 5.11 Converged solution for projectile base flow problem using the TVD scheme 0.2 0.0 -0.2 -0.4 -0.6 0.0 -0. -0.4 0.37 CHAPTER VI PRELIMINARY TESTS AND MODIFICATIONS VI. 1 The General Procedure Prior to coupling the Adaptive Grid code with the flow solvers to solve the projectile flow problem, it was found beneficial to modify the adaptive grid generation method. The following sections include an example showing the control over orthogonality and then describes two modifications made to the adaptive grid scheme to improve the general technique. It should be noted that each of the modifications is made to improve the method for the present flow problem and geometry and may not be necessary or useful in other problems. For each example the initial grid network used was that of Figure 5.5, or that of 5.10 if the base flow configuration is used in the example. In the adaptive grid code, the control function P was set equal to a constant so that the clustered points along each C coordinate line in the initial grid would tend to spread and equal spacing would result. The control function Q was defined using eqs. (4.3.2) with the distribution of eq. (4.3.3) so that the same point distribution along the r7 coordinates in the initial grid would be retained in the grid obtained using the Adaptive Grid code. Each of the figures used in the examples show only the points in the reduced grid network. 62 VI.2 Orthocronality In order to demonstrate the effect of enforcing orthogonality, two cases were considered, one with the parameter A = 0.0, the second with = 0.5. For the present projectile, the most difficult area to obtain an orthogonal grid is over the projectile nose. This is due to the sharp nose configuration which results in the two boundaries, the projectile surface and the upstream axis of symmetry, meeting at an obtuse angle. The two sets of coordinates, which run parallel to the boundaries, then tend to intersect at similar angles. This result is shown in Figure 6. la, which shows an enlarged view of the grid in the nose section that resulted when A = 0.0. To improve the orthogonality in this region, another case for A = 0.5 was run, and as the results of Figure 6. 1b indicate, the orthogonality of the grid coordinates near the nose was in general increased. For all the solutions to the transonic flow problem using the adaptive grid technique reported in Chapter VIII, A was set equal to 0.5. VI.3 Curvature One inherent result of grids obtained as solutions to equations based on the Laplacian operators is the effect of curved boundaries on the spacing of interior coordinate lines. In general, boundaries that are convex will attract coordinate lines and those that are concave will repel coordinate lines. In the context of the adaptive grid generation scheme presented here, this effect will alter the grid spacing dictated by the control functions, either increasing or decreasing the spacing that would result in the absence of curvature effects. The severity of the boundary curvature's influence on the grid is best seen in solving Figure 6.la Grid network using X=0.0 Figure 6.1b Grid network using X=0.5 64 the adaptive grid generation equations for the base flow grid network. Recall that thel control function Q is defined so that the original grid point spacing in the normal direction would be retained. However, the point spacing resulting as a solution to the adaptive grid equations differs along each 7 coordinate line as shown in Figure 6.2a. The attraction of points is most prominent near the sharp corner at the projectile base where the effect is large due to the high curvature at the corner. A simple and effective technique to eliminate the effect of curvature can be developed by considering the solution of Laplace's equations between two concentric circles C1 and C2 with radii R1 and R2 respectively. V2 = 0 (6.2.1) V2 = 0 Solving Laplace's equation is equivalent to solving the adaptive grid equations with constant control functions and neglecting orthogonality (i.e. A = 0.0). In using constant control functions, we expect the coordinates to be equally spaced in each direction. It is advantageous to write the equations in polar coordinates (r, ), and it is sufficient here to seek a solution in which q is a function of 0 only and 7 is a function of r only. The coordinate lines then align themselves with the 0 coordinate direction and the ri coordinate lines become aligned with the radial coordinate r. The problem can then be formulated as Figure 6.2a Grid network obtained without corrections for curvature Figure 6.2b Grid network obtained with corrections for curvature -0 d02 d2 i 1 do + - =0 dr2 r dr (6.2.2) = 0 on C1 and C2 S= R1 on C1 = R2 on C2 The boundary conditions corresponds to uniform spacing along the coordinate lines coincident with the bounding concentric circles. The spacing along each coordinate direction is the inverse of the first derivatives O and Tr which by integrating eqs. (6.2.2) are found to be 8 = a b (6.2.3) Tr r where a and b are constants of integration. These results indicate that the spacing in the direction is uniform while the spacing in the ?I direction is smaller near the inner circle C1 and larger near the outer circle C2 The term in the governing equation responsible for this result, designated here as K, is 1 dT7 K = - (6.2.4) r dr which can be viewed as the curvature of the intersecting coordinate line, (1/r), multiplied by the inverse of the spacing along the tI 67 coordinate line. The effects of the curved coordinates on the spacing along the r coordinate lines can be eliminated for this problem simply by subtracting K frame the governing equation for Y7, V2 = 0 (6.2.5) S1 d? V2 =0 r dr This technique can be extended to the general adaptive grid equations by developing the curvature term K for each curvilinear coordinate line in the computational domain and then subtracting them fram eqs. (2.2.7). The term needed to cancel the effects of a curved rP coordinate line on the spacing along a coordinate line is comprised of the curvature of the i) coordinate line and the spacing along the C coordinate line. The curvature p of a n coordinate line is defined as d p = I (6.2.6) ds where r is a unit tangent to the rI coordinate and is in the computational plane -_ ) (x,,y) --Z -(6.2.7) Jaz where a is defined in eq. (2.3.5). By using the two relationships S = / d_ i d (6.2.8) ds s7 dt7 the curvature p can be shown to be 68 =( yx -- xLI y-17t (6.2.9) a3/2 By multiplying this result with the inverse of the spacing along a coordinate (1/7-y), the following expression is obtained Sy~x x~y,, K1 =f? a? 17 (6.2.10) j7 ct3/ 2 which is the term that is subtracted from the adaptive grid equation for in order to cancel the effects of curved t7 coordinate lines. An expression for the influence of curved coordinates on the spacing along P7 coordinate lines can be obtained in a similar manner and is y~x^ x~y^ K2 = (6.2.11) j/a 'Y- 3/2 Thompson et al. [15] have derived similar expression, but do not implement them in the same way. It should be noted that in the computational plane, the terms K1 and K2 contain second derivatives of the variables x and y, and when eqs. (2.2.7) are transformed onto the computational plane, the second order central difference approximations to these terms depend explicitly on the point (xi j, Yi,j). However, in implementing the Newton-Raphson iterative solution procedure, the terms K1 and K2 are considered part of the source terms. Thus, the right hand side of eqs. (2.3.7) becomes UJ P S1 K1J3 -'- P Klj (6.2.12) aJ Q r, S2 K2J3 This is equivalent to neglecting the dependence of the finite difference approximation on the point (xi,j, Yi,j) in obtaining eqs. (3.1.5) for the Newton-Raphson iterative procedure. The J3 term in eqs. (6.2.12) appears in the transformation of eqs. (2.2.7) into the ccaputational plane. The adaptive grid generation equations, modified in this manner, will produce the grid shown in Figure 6.2b. In contrast to Figure 6.2a, it is clear that the effects of the curved boundaries have been successfully eliminated. VI. 4 Internal Grid Boundaries One other problem that arose in the preliminary testing of the adaptive grid generation technique is the rapid upstream bending of the ri coordinate lines emanating from the secant portion of the projectile surface. Figure 6.3a shows a grid network obtained using the adaptive grid generation equations with the curvature terms included and with x = 0.5. As can be seen, the Pi coordinate lines emanating frcm the secant portion of the projectile quickly bend forward to fill the upstream region of the domain which results in a very skew grid. This characteristic of the grid does not appear as a problem until a solution using the Beam and Warming scheme is sought on such a grid. Figure 6.3b shows a time history of the solution obtained on such a grid for Mach 0.96 and as can be seen the pressure coefficient oscillates in time over Figure 6.3a Grid network upstream of the projectile MACH 0.91 0.2 0.0 -0.2 -0.4 -0.6 I o experiment C- computed Figure 6.3b Lack of convergence over secant ogive with Beam and Warming scheme 71 the secant portion of the projectile. The solution did not converge and instead the oscillations continued and eventually effected the solution downstream. In order to prevent the upstream bending of the ri coordinates, one coordinate line in the initial grid network is designated as an internal boundary. This rt coordinate line remains fixed in space during the iterativesolution procedure, thus no other rl coordinate lines will pass through it. However, the points defining the internal boundary coordinate line are allowed to move along the coordinate so that they remain consistent with the spacing of the points along the adjacent r, coordinate lines. This is accomplished by updating the arclength distribution along the internal boundary after each iterative sweep in which the new distribution is defined as the average of the arclength distributions along the two adjacent rY coordinate lines. Figure 6.3c shows the grid obtained when the 12th q coordinate line from the projectile nose is designated as an internal boundary and as can be seen it prevents the coordinate lines from bending upstream. It has been chosen to lie in a region in which no major grid clustering is expected so that it will not effect the adaptation of the grid network. The length and number of points used along the C coordinate lines in each section may result in an abrupt change in the point spacing along the coordinate lines passing through the boundary as can be seen in Figure 6.3c. Therefore the control function P, controlling the spacing along the coordinate lines is locally modified near the internal boundary so that the grid point spacing varied smoothly through the internal boundary. This can be accomplished by requiring the normalized rate of change of grid point spacing along a coordinate line (sd/sC) Figure 6.3c Grid network with internal boundaries Figure 6.4 Grid points near internal boundary 73 be identical at the two points adjacent to the internal boundary. Referring to Figure 6.4, if the h coordinate line designated as an internal boundary contains the ith point along a coordinate line, we then require ( --] =s-- ] (6.3.1) s I i-i= s" I i+1 To implement this requirement into the solution procedure, consider the one dimensional equation for adaptation along a C coordinate, eq. (2.3.8) which is written here as P s PC __ (6.3.2) P s The left hand side of eq. (6.3.2) appears in the source terms of the two dimensional equations for adaptation, eqs. (2.3.7). Substitution of eq. (6.3.2) into the source terms for a grid point just ahead of the internal boundary results in the modified source terms at the point i-1 f s 1 633 (Si)i-l.'j = ai-ljJi-lj I (6.3.3) The modified source terms for the point just behind the internal boundary are obtained by replacing i-1 with i+1 in eq. (6.3.3). In order to apply the iterative solution technique a value for the terms (sc/sc) in eqs. (6.3.3) for each point adjacent to the internal boundary will be needed. Their final values are not known a priori since they will depend on the final position of the points in each section, therefore they are set to the current value [s 1_ (si-2 -2si + si+2) (6.3.4) s Ji-l,j (si+2-si-l)/2 The grid that is obtained when using this technique is shown in Figure 6.5 in which it can be seen that the abrupt change in spacing across the internal boundary of the grid in Figure 5.3c has been reduced. The use of internal boundaries provides a simple and effective way to control the grid points in the domain. In fact, it was found helpful to also designate the Y coordinate line which separates the region over the sting from that over the projectile as an internal boundary to keep the 7 coordinate lines that originally began on the projectile surface from moving over the sting. Figure 6.5 Smoothed grid spacing near the internal boundary CHAPTER VII THE SELF-ADAPrIVE CaUATIONAL METHOD The adaptive grid generation technique can be incorporated naturally into the solution process of the flow solvers since the flow solver algorithms are time-marching. After a specified number of time steps, the grid network can be adapted to the control functions based on the current values of the flow variables. Thus the time-marching algorithm of the two flow solvers used here, the Thin Layer and TVD codes, proceeds as usual with intermittent calls to the Adaptive Grid code to update the grid network. In this approach, however, the relocation of the grid points in the physical domain due to the grid adaptation must be considered. The values of the flow variables at any time step are defined at the current grid point locations and any movement of the points will therefore distort the solution. This effect is critical for unsteady flow problems since the grid point movement will effect the solution accuracy. For steady flow problems, the grid point movement may effect the convergence rate and if it is large enough, it may cause the solution to diverge. There are two basic approaches to account for the grid point movement resulting when adapting the grid network. The first is to view the grid network as a moving grid and include the time metrics in the coordinate transformation of the Navier-Stokes equations. These metrics can be written in the computational plane as 77 Ct = CxXr + y Yr (7.1) I?t = 7xXr + 'yYr The quantities x, and y, can be approximated by using a backward difference expression x- = (Xn+l-xn)/At (7.2) Yr. (yn+l-yn)/At (7.2) two in which xn and yn are the grid point locations of the previous grid network and xn+1 are yn+l are the values of the adapted grid network. This approach requires the grid to be adapted every time step which is inefficient for steady flow problems. Furthermore, it has been shown that special procedures must be used in evaluating the time derivative of the Jacobian [20] in order to maintain the conservative property of the finite difference schemes. Another approach, used here, is to interpolate the flow variables from the previous grid network onto the adapted grid network. This approach can be applied to both steady and unsteady flow problems and does not require the grid to be adapted every time step. Thus this approach can be used for both steady and unsteady flow problems. For the two- dimensional grid networks considered here, a linear interpolation is currently used. Each grid cell in the previous grid network is divided into two triangles. Once a grid point in the adapted grid network is located within a triangle of the previous grid network, the value of the flow variables at that point can be determined by interpolation. Referring to Figure 7.1, the value of a flow variable q at the point in the adapted grid network can be determined from the values of the flow / A1 . A, . 0 Previous grid network * current grid network Figure 7.1 Flow variable interpolation 79 variables at the vertices of the triangle as (q1A1 + q2A2 + q3A3) q = (7.3) (A1 + A2 + A3) where A,, A2 and A3 are the areas of the triangles formed with the grid point in the adapted grid network. The basic algorithm for incorporating the adaptive grid generation technique into both the Thin Layer code and the TVD code consists of the following steps: 1. Input an initial grid network into the flow solver code and begin the time-marching algorithm. 2. Advance the solution a specified number of time steps, and submit the current grid network to the Adaptive Grid code. 3. Obtain an adapted grid network using the Adaptive Grid code. 4. Interpolate the flow variables onto the adapted grid network which now becomes the current grid network. 5. Repeat steps 2 through 5 until the solution converges for steady flow problems or repeat the steps until a specified time is reached for unsteady flow problems. This process which couples the flow solver and the adaptive grid generation technique is referred to here as the self-adaptive computational technique. For the two problems considered here, the axisymmnetric transonic flow over a projectile with sting and the transonic projectile base flow problem, only steady flow solutions are sought. Since the adapted grid network depends on the flow variables via the control functions, the convergence of the solution implies also that the grid network ceases to change with time. In all the examples of the following chapters the grid networks of Figures 5.4 and 5.10 are used as 80 the initial grid. For the purpose of adapting the grid network, the source terms corrected to eliminate curvature effects, eqs. (6.2.12) were used in the adaptive grid generation equations. In solving the projectile problem with sting two internal t7 coordinates were designated as internal boundaries. The first is over the secant ogive as shown in Figure 6.5, the second separates the grid network over the projectile from that over the sting and is located at approximately x=7.0. In solving the projectile base flow problem, no internal boundaries were used. CHAPrER VIII RESULTS AND DISCUSSION VIII. 1 Choices for the Control Functions In applying the self-adaptive computational technique it is necessary to define the control functions for adaptation. For the transonic projectile problem considered here both shocks and expansion waves will occur approximately normal to the projectile surface and require refined grid point spacing in the coordinate direction for adequate resolution. The original choice for the function f in the control function P to obtain this clustering was the pressure gradient. However, it was found in the numerical experiments that the expansion waves occurring at the secant ogive/cylirnder and cylinder/boattail junctures required smaller spacing than the shocks for good resolution. The pressure gradient was largest in the shocks and thus smaller spacing occurred in the shock regions rather at the expansion waves. A better choice for the control function was found to be the curvature of the pressure distribution along the streanwise coordinate direction, as2 f = ----- I 8S I(8.1.1) ( 1 + (an )2)3/2 as 82 where the derivatives are with respect to the arclength s along each coordinate line. This function is largest at the expansion waves where the pressure gradient changed sign. At the base and peak of a shock eq. (8.1.1) also increased but it is not as large as the value near the expansion waves. Thus the smallest spacing will occur near the expansions and yet the grid points will also cluster, to a lesser extent, in the vicinity of shocks. In the normal direction, the important area requiring increased resolution is the boundary layer and the logical choice for the function g in the control function Q is the velocity gradient. However, it was found that the control function based on the velocity gradient resulted in too many points being placed in the boundary layer and led thus to a rapid stretching of the grid point spacing along the Y7 coordinate lines. Another choice for the function g, the exponential of the velocity gradient g = exp as (8.1.2) where u is the velocity magnitude and s measures along the r, coordinate line, was found to yield better results. Figure 8.1 shows a plot of the control function Q along a typical S,, coordinate line. The control function resulting from the exponential of the velocity gradient is similar to the control function corresponding to the grid point distribution in the fixed grid networks which is known to give good results provided enough points are used. The control function obtained using the velocity gradient however is large for many grid points and, o velocity gradient * prescribed distribution o exponential of the velocity gradient 20 40 77 60 Figure 8.1 Comparison of control functions MACH 0.96 o experiment - computed Figure 8.2 Lack of convergence using a self-adaptive computational technique with Beam and Warming scheme Q (1o) 3.0 1.5 0.2 0.0 -0.2 -0.4 -0.6 84 thus, results in too many points being clustered into the boundary layer region at the expense of adequate resolution in other areas. The minimum and maximum grid point spacing along the t coordinate lines is 0.00002 and 6.0 respectively. In this direction it is important to have the small spacing at each point along the projectile surface in order to adequately resolve the viscous sublayer region. The parameter a of eq. (4.1.4) was determined for each j coordinate line so that the small spacing would occur for each point on the projectile surface. However, in the coordinate direction, the shocks and expansion waves are not expected to reach the outer boundary and the small spacing needed along the projectile surface is not required near that boundary. With the minimum and maximum grid point spacing prescribed as 0.02 and 0.22, respectively, a is determined only for the coordinate line coincident with the projectile surface and that value is used for all the remaining coordinate lines. This had the effect of increasing the minimum spacing as the magnitude of the shocks and expansion waves diminished near the outer boundary. VIII. 2 Results Using the Self-Adaptive Computational Technique: The Beam and Warming Scheme The first application of the self-adaptive coputational technique is with the Beam and Warming scheme for the calculation of the projectile flow with sting. The Mach number is 0.96 and the Reynolds number is 76,000, a case for which experimental data is available. The time step used in the calculation is 0.1 and the explicit and implicit dissipation terms are 2 and 4 times the time step respectively. The initial grid network used is that of Figure 5.5 with IM=70 and JM=60. In the first 85 attempts to use the adaptive grid technique, the grid network was adapted every 200 time steps. A series of calculated pressure coefficients obtained with this approach are shown in Figure 8.2. As can be seen, the pressure coefficient oscillates widely over the projectile surface and shows no indication of converging. In a series of attempts to obtain better results, the time step was decreased to 0.05, the dissipation terms were doubled and the number of time steps between adaptation were both increased to 500 and reduced to 50. However, similar results were obtained in each case. In the next case tried, the number of points in the streamwise direction was increased to 90 points, thus IM=90, and the grid network was adapted every time step. The results for this case are nuch better, however a converged solution was not obtained. Figure 8.3a shows the computed pressure coefficient at the times indicated, and as can be seen, the solution appears to have converged everywhere except in the shock boundary layer interaction regions. The shocks propagate upstream and downstream in a cyclic fashion. Figure 8.3b shows the same phenomenon occurring at a later time. In fact, the solution was continued for 8000 time steps and the propagation of the two shocks showed no indication of diminishing. Figure 8.4a shows the reduced grid network which resulted at 5000 time steps. As is evident, the grid point spacing in the q coordinate direction is clustered near the projectile surface with a distribution similar to that of the hyperbolic grid in Figure 5.4. The spacing varies somewhat in the shock and expansion wave regions, most likely as a result of enforcing orthogonality. The clustering in the streamwise direction MACH 0.96 5a a o Us wes- .^ ite, o experiment computed Figure 8.3a Computed solution between 4000 and 5000 time steps with Beam and Warming scheme Cp MACH 0.96 0.2 h 6000 -0.4 o experiment -0.6 computed Figure 8.3b Computed solution between 5000 and 6000 time steps with Beam and Warming scheme Cp 0.2 0.0 -0.4 -0.6 -0.2 Figure 8.4a Adapted grid network at 5000 time steps (90 x 60) Figure 8.4b Adapted grid network at 6000 time steps (90 x 60) 88 clearly indicates a response to the control function P, which is based on the curvature of the pressure function. Clustering occurs at both junctures on the projectile surface and the adaptation to the shock structures on the cylinder and boattail are also quite evident. Note that the smallest spacing appears to be located in the shocks well above the projectile surface. This occurs because there is more competition for the small grid spacing on the projectile surface since both the shocks and expansion waves are strong in that region. Figure 8.4b shows the grid network corresponding to a later time. In this Figure and those following, the coordinate lines are not shown in order to gain a better view of the grid spacing near the projectile surface. The grid configuration is basically the same as that of Figure 8.4a, however, the position of the shocks as indicated by the adaptation has changed. The shock on the cylinder has propagated downstream and the boattail shock has moved upstream. Referring to the shock positions evident in Figures 8.3a and 8.3b it can be seen the clustering in the grid network is successfully following the shock patterns. In an attempt to improve these results the time step was reduced to 0.05 and the dissipation terms were again increased, however, similar results were obtained. Also the convergence criterion 6 of eq. (3.1.12) used in the iterative solution procedure for the grid generation equation was varied. At the value of 6=0.005, only one iteration was required for the grid network to converge. The rapid convergence occurred because the grid was updated every time step and the solution, and consequently the control functions, did not change much between grid adaptations. When the criterion was reduced to 6=0.0025, 1 to 3 iterations were required for 89 the grid to converge at each time step, however it had no effect on the cyclic propagation of the shocks. Although a converged solution was not obtained using the adaptive grid generation technique with the Beam and Warming scheme, there are still some encouraging results with respect to the adaptive grid generation procedure. The maximum and minimum grid point spacing occurring along the projectile surface in the adapted grid network were within 5 percent of their prescribed values. The average minimm spacing along all the rj coordinate lines was within 3 percent of the prescribed value of 0.00002 with a maximum deviation of 10 percent. The adapted grid networks shown in Figures 8.4a and 8.4b show that the adapted grid generation equations are capable of responding reliably and predictably to the chosen control functions and yield the prescribed maximum and minimum grid point spacing. The adaptive grid generation equations, thus, have been shown to provide a well adapted grid network for the transonic flow problems provided that a good choice is made for the control functions. However, the convergence of the solution was effected by the use of the adaptive grid generation procedure with the Beam and Warming scheme. One reason for the failure of the solution to converge may be due to interpolation. Error introduced due to interpolating the flow variables onto each successive adapted grid network may continuously perturb the solution, causing it to oscillate. Another cause of this phenomenon may be the form of the dissipation terms in the Beam and Warming scheme. The added dissipation terms contain derivatives with respect to the ccuputational coordinates and were added to the scheme after the governing equations were transformed. Since these terms are a function of the transformed 90 coordinates any change in the grid network used will change the dissipation terms and consequently alter slightly the transformed governing equations. Thus the continuous adaptation of the grid network continuously alters the dissipation terms and may be the reason the solution does not converge. VIII. 3 Results Using the Self-Adaptive Cocmputational Technique: The TVD Scheme Since the solutions obtained using the Beam and Wanrming scheme did not converge, another scheme, the TVD scheme was also used. Recall that the primary difference between the two schemes is the dissipation terms added to the governing equations. The initial grid network used is that of Figure 5.4 with IM=70 and JM=60. The time step used is 0.1 and the grid was adapted every 200 time steps. For Mach 0.96, the solution converged easily within 2000 time steps. Only one iteration was needed for the grid network to converge at 1800 and 2000 time steps which also indicates a converged solution. Figure 8.5a shows the computed pressure coefficient using the adapted grid generation procedure compared to the solution obtained on the fixed grid network of Figure 5.4. The solutions agree well with experimental data over most of the projectile surface, except for the predicted position of the shock over the cylinder. The position predicted using the adaptive grid generation technique is centered on the cylinder whereas the position predicted using the fixed grid network is slightly upstream of the center. Figure 8.5b shows a shadowgraph for corresponding flow conditions in which it can be seen that the experimentally predicted shock location is centered over the cylinder portion of the projectile. Thus, the use of the adaptive grid generation technique improved the solution accuracy. Figure 8.5c shows 91 the adapted grid network corresponding to the converged solution. Again, clustering near the secant ogive/cylinder and cylinder/boattail junctures is evident as well as the adaptation to the shocks. In comparison to the fixed grid network of Figure 5.4 it is evident that the grid spacing along the projectile surface is refined at the center of the cylinder portion due to the adaptation and is responsible for the increased solution accuracy. Two other cases were run using the identical procedure used for the case of Mach 0.96. Figure 8.6a shows a comparison of the solution obtained using both the fixed grid network and the adaptive grid generation technique for Mach 0.91, another case for which experimental data is available. The only improvement obtained using the adaptive grid generation procedure is an increase in the peak of the shock occurring over the cylinder portion of the projectile. However, the peak still does not reach the experimentally predicted value. In an attempt to increase the accuracy in that region, 10 more points were added in the coordinate direction, thus IM=80, however no further improvement in the solution obtained with the adaptive grid generation technique was evident. The results thus may be the best that the TVD scheme can yield given the current grid configuration and the shortened pressure peak is probably due to the larger dissipative error associated with the TVD scheme. Figure 8.6b shows the grid network corresponding to the converged solution. Again grid point clustering in response to the curvature of the pressure is evident, and as can be seen by the adaptation of the grid network, both shocks occur directly downstream of the expansion waves. MACH 0.96 s2- C' / ,/, I K7 0.4- - adopted grid . fixed grid Figure 8.5a Comparison of computed solutions on fixed and adaptive grid networks with TVD scheme T j -B k.. .. Figure 8.5b SOCBT shadowgraph for Mach 0.96 (reprinted from reference 24) 0.2 0.0 -0.2 Figure 8.5c Adapted grid network for converged solution for Mach 0.96 Cp MACH 0.91 0.2[ o 0.2 adapted grid .-- fixed grid -0.6 Figure 8.6a Comparison of computed solutions on fixed and adapted grid networks with TVD scheme |

Full Text |

16
spacing along a coordinate line. If the control functions change rapidly along a coordinate, the grid spacing will consequently change rapidly. This smoothing thus can be governed through the definition of the control functions. Another type of smoothness, which is inherent in the elliptic equations, is that between the grid point spacing of adjacent coordinates. If the grid point distribution along two adjacent coordinates varies, the elliptic equations will tend to dampen the variation. Equations (2.2.7) govern the grid coordinates Â£ and rÂ¡ and will be elliptic as long as the parameter A is not too large. They constitute a boundary value problem and it is therefore necessary to define the coordinates along the boundaries. It is also necessary to adapt the boundary points in a consistent manner with the interior points so that the boundary points remain aligned with the interior points. For this purpose, a one dimensional functional analogous to equation (2.2.1) or (2.2.2) is defined to measure adaptation along a boundary coordinate, (2.2.8) where s is the arclength along the boundary coordinate. There is no consideration for orthogonality and scaling is unnecessary. Applying the one dimensional Euler-Lagrange equation yields a second order differential equation, Â£ss Â£ s Ps/P0 (2.2.9) Note that a similar equation can be obtained for boundary adaptation along a rÂ¡ coordinate by replacing f in eq. (2.2.9) with r? and P with Q. This one dimensional equation, together with eqs. (2.2.7), are the 88 clearly indicates a response to the control function P, which is based on the curvature of the pressure function. Clustering occurs at both junctures on the projectile surface and the adaptation to the shock structures on the cylinder and boattail are also quite evident. Note that the smallest spacing appears to be located in the shocks well above the projectile surface. This occurs because there is more competition for the small grid spacing on the projectile surface since both the shocks and expansion waves are strong in that region. Figure 8.4b shows the grid network corresponding to a later time. In this Figure and those following, the Â£ coordinate lines are not shown in order to gain a better view of the grid spacing near the projectile surface. The grid configuration is basically the same as that of Figure 8.4a, however, the position of the shocks as indicated by the adaptation has changed. Ihe shock on the cylinder has propagated downstream and the boattail shock has moved upstream. Referring to the shock positions evident in Figures 8.3a and 8.3b it can be seen the clustering in the grid network is successfully following the shock patterns. In an attempt to improve these results the time step was reduced to 0.05 and the dissipation terms were again increased, however, similar results were obtained. Also the convergence criterion 5 of eq. (3.1.12) used in the iterative solution procedure for the grid generation equation was varied. At the value of 5=0.005, only one iteration was required for the grid network to converge. The rapid convergence occurred because the grid was updated every time step and the solution, and consequently the control functions, did not change much between grid adaptations. When the criterion was reduced to 6=0.0025, 1 to 3 iterations were required for Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy A SELF-ADAPTIVE COMPUTATIONAL METHOD FOR TRANSONIC TURBULENT PROJECTILE AERODYNAMICS By Christopher William Reed August, 1987 Chairman: Chen-Chi Hsu Major Department: Engineering Sciences An adaptive grid generation procedure is developed for viscous flow problems. The equations governing the adaptation are based on a variational statement resulting in a set of elliptic equations in which adaptation can occur independently in each coordinate direction. The method allows for explicit control of adaptation and orthogonality while grid smoothness is implicit in the elliptic equations. The adaptive grid generation equations retain a simple relationship between the control functions and the grid point spacing, are capable of efficiently providing the extremely refined mesh in the boundary layer regions and can maintain a smooth grid network in complex geometries. A self-adaptive computational technique has been developed by coupling the adaptive grid generation procedure with a thin layer Navier-Stokes code and a TVD code. In solving an axisymmetric turbulent vi CHAPTER VIII RESULTS AND DISCUSSION VIII.1 Choices for the Control Functions In applying the self-adaptive computational technique it is necessary to define the control functions for adaptation. For the transonic projectile problem considered here both shocks and expansion waves will occur approximately normal to the projectile surface and require refined grid point spacing in the Â£ coordinate direction for adequate resolution. The original choice for the function f in the control function P to obtain this clustering was the pressure gradient. However, it was found in the numerical experiments that the expansion waves occurring at the secant ogive/cylinder and cylinder/boattail junctures required smaller spacing than the shocks for good resolution. The pressure gradient was largest in the shocks and thus smaller spacing occurred in the shock regions rather at the expansion waves. A better choice for the control function was found to be the curvature of the pressure distribution along the streamwise coordinate direction, f = & 2 ( 1 + (** )2)3/2 as (8.1.1) 81 62 VI.2 Orthogonality In order to demonstrate the effect of enforcing orthogonality, two cases were considered, one with the parameter A = 0.0, the second with A = 0.5. For the present projectile, the most difficult area to obtain an orthogonal grid is over the projectile nose. This is due to the sharp nose configuration which results in the two boundaries, the projectile surface and the upstream axis of symmetry, meeting at an obtuse angle. The two sets of coordinates, which run parallel to the boundaries, then tend to intersect at similar angles. This result is shown in Figure 6.1a, which shows an enlarged view of the grid in the nose section that resulted when A = 0.0. To improve the orthogonality in this region, another case for A = 0.5 was run, and as the results of Figure 6.1b indicate, the orthogonality of the grid coordinates near the nose was in general increased. For all the solutions to the transonic flow problem using the adaptive grid technique reported in Chapter VIII, A was set equal to 0.5. VI.3 Curvature One inherent result of grids obtained as solutions to equations based on the Laplacian operators is the effect of curved boundaries on the spacing of interior coordinate lines. In general, boundaries that are convex will attract coordinate lines and those that are concave will repel coordinate lines. In the context of the adaptive grid generation scheme presented here, this effect will alter the grid spacing dictated by the control functions, either increasing or decreasing the spacing that would result in the absence of curvature effects. The severity of the boundary curvature's influence on the grid is best seen in solving 69 aJ P Si = ^ KiJ3 (6.2.12) s2 = qJ Q n Q k2J3 This is equivalent to neglecting the dependence of the finite difference approximation on the point (xj_f j, yifj) in obtaining eqs. (3.1.5) for the Newton-Raphson iterative procedure. The J3 term in eqs. (6.2.12) appears in the transformation of eqs. (2.2.7) into the computational plane. The adaptive grid generation equations, modified in this manner, will produce the grid shown in Figure 6.2b. In contrast to Figure 6.2a, it is clear that the effects of the curved boundaries have been successfully eliminated. VI.4 Internal Grid Boundaries One other problem that arose in the preliminary testing of the adaptive grid generation technique is the rapid upstream bending of the rÂ¡ coordinate lines emanating from the secant portion of the projectile surface. Figure 6.3a shows a grid network obtained using the adaptive grid generation equations with the curvature terms included and with A = 0.5. As can be seen, the rÂ¡ coordinate lines emanating from the secant portion of the projectile quickly bend forward to fill the upstream region of the domain which results in a very skew grid. This characteristic of the grid does not appear as a problem until a solution using the Beam and Warming scheme is sought on such a grid. Figure 6.3b shows a time history of the solution obtained on such a grid for Mach 0.96 and as can be seen the pressure coefficient oscillates in time over 15 Euler-Lagrange equations. This violates the variational principal, but it has been shown [13] that the use of local scaling will produce better grids when there are large changes in the grid spacing. After substituting in the scales, the total functional Ip becomes It= [ { Pl V^:VÂ£ + QlV>? Vr? + ajL (v^.Vr?)2 }dxdy (2.2.6) After applying the Euler-Lagrange equations for Â£ and rÂ¡ to the total functional I c + c + yp-.ve xx + + ^ (*?xÂ£xx+2r7x7yÂ£xy+riyyy + + (v xÂ£ y+r) yÂ£ x) 7 xy+ (Â£ x7 x+2Â£y7y) 7yy) (2.2.7) 7XX + 7yy + + (?x7x+2Â£y7y) Â£xy + * ( (2Â£x7x+Â£y7y) ?xx+ (rlx^y+riy^x) Â£xy + 2 p ^x7xx+2^x^y7xy + Cy7yy)=0 The subscript L has been dropped from the local scales in the differential equations since they are no longer needed to indicate their assumed invariance. The Â£ and rÂ¡ coordinates that satisfy these equations will yield a grid with the desired properties defined in the functionals of eqs. (2.2.1), (2.2.2) and (2.2.3). It should be noted that there was no explicit functional to maintain grid smoothness. There are two types of smoothness that should be differentiated. One type, which can be controlled by the control functions, is smoothness of the grid point 89 the grid to converge at each time step, however it had no effect on the cyclic propagation of the shocks. Although a converged solution was not obtained using the adaptive grid generation technique with the Beam and Warming scheme, there are still some encouraging results with respect to the adaptive grid generation procedure. The maximum and minimum grid point spacing occurring along the projectile surface in the adapted grid network were within 5 percent of their prescribed values. The average minimum spacing along all the r? coordinate lines was within 3 percent of the prescribed value of 0.00002 with a maximum deviation of 10 percent. The adapted grid networks shown in Figures 8.4a and 8.4b show that the adapted grid generation equations are capable of responding reliably and predictably to the chosen control functions and yield the prescribed maximum and minimum grid point spacing. The adaptive grid generation equations, thus, have been shown to provide a well adapted grid network for the transonic flow problems provided that a good choice is made for the control functions. However, the convergence of the solution was effected by the use of the adaptive grid generation procedure with the Beam and Warming scheme. One reason for the failure of the solution to converge may be due to interpolation. Error introduced due to interpolating the flow variables onto each successive adapted grid network may continuously perturb the solution, causing it to oscillate. Another cause of this phenomenon may be the form of the dissipation terms in the Beam and Warming scheme. The added dissipation terms contain derivatives with respect to the computational coordinates and were added to the scheme after the governing equations were transformed. Since these terms are a function of the transformed Figure 3.1 High aspect ratio cell Figure 3.2 Initial grid network 19 = *r,2 + Yb2 7 = XÂ£2 + y^2 P = + xrjYr, A' = A/J and Si = P^Ja P PnJ/9 s2 = Q Qt7J7 + Q (2.3.6) It is important to note that in inverting the equations the derivative of P in the r? direction appears even though P was intended to control the spacing in the Â£ direction. This occurs because the contravariant base vector VÂ£ actually indicates the spacing in the direction normal to lines of constant Â£ rather than in the Â£ direction. If the grid is orthogonal, the two directions, VrÂ¡ and VÂ£, are equivalent and the r? derivative will not have any influence since p will vanish. It has been found in numerical experiments that eliminating this term from the equations has a negligible effect and consequently it is dropped from the equations. The source terms S1 and S2 in eqs. (2.3.5) then become P^Ja Qt7J7 (2.3.7) S2 = Q REFERENCES 1. Kutler, P, "A Perspective of Theoretical and Applied Computational Fluid Dynamics," AIAA Journal, Vol. 23, NO. 3, pp. 328-341, 1984. 2. Thompson, J.F., Thames, F.C., and Mastin, C.W., "Automated Numerical Generation of Body-Fitted Curvilinear Coordinate System for Field Containing Any Number of Arbitrary Two-Dimensional Bodies," J. Comp. Phys., Vol. 15, pp. 299-319, 1974. 3. Mastin, C.W., "Error Induced By Coordinate Systems," Numerical Grid Generation. Ed. J.F. Thompson, North-Holiand, p. 41, 1981. 4. Hsu, C.C., and Chang, T.H., "On a Numerical Solution of Incompressible Turbulent Boundary Layer Flow," Finite Element Flow Analysis. Ed. T. Kawai, Univ. of Tokyo Press, Tokyo, pp 219-226, 1982. 5. Saltzman, J. and Brackbill, J.U., "Application and Generalization of Variational Methods for Generating Adaptive Meshes," Numerical Grid Generation. Ed. J.F. Thompson, North-Holland, pp. 865-884, 1981. 6. Hsu, C.C. and Tu, C.G., "An Adaptive Grid Generation Technique Based on Variational Principals for Transonic Projectile Aerodynamic Calculation," To Appear in Int. J. Numerical Methods in Fluids. 7. Nakahashi, K., and Deiwert, G.S., "A Self-Adaptive Grid Method with Applications to Airfoil FLow," AIAA Paper 85-1525. 8. Dywer, H.A., Smooke, M.D. and Kee, R.J., "Adaptive Gridding for Finite Difference Solutions to Heat and Mass Transfer Problems," Numerical Grid Generation. Ed. J.F. Thompson, North-Holland, New York, p. 339, 1982. 9. Gnoffo, P. A., "A Vectorized, Finite Volume, Adaptive-Grid algorithm for Navier-Stokes Calculations," Numerical Grid Generation. Ed. J. F. Thompson, North-Holland, New York, p. 819, 1981. 10. Thompson, J.F., "A Survey of Dynamically-Adaptive Grids In the Numerical Solution of Partial Differential Equations," Applied Numerical Mathematics 1. North-Holland, New York, 1985. 11. Rai, M.M. and Anderson, D.A., "Grid Evolution in Time Asymptotic Problems," J. of Comp. Rhys., Vol 43, pp. 327-344, 1981. 105 6 an adaptive grid generation technique. First is the definition of the control functions, which are the part of the grid generation equations that dictate the need for adaptation. They are dependent upon the solution and thus provide the link between the grid generation and the governing equations. Also, a method is needed to account for the effect of the moving grid on the solution. Since the solution is defined at each point, any movement will alter the solution. The proposed adaptive grid generation technique is used to adapt the grid for a transonic axisymmetric projectile flew problem which is solved using both an implicit factorized algorithm and a TVD scheme for the thin layer Navier Stokes equations. The constraint of axisymmetric flow reduces the domain to two independent spatial variables, thus the grid generation equations are developed in two dimensions. The technique, however, can be readily extended to three dimensions. 87 Figure 8.4a Adapted grid network at 5000 time steps (90 x 60) Figure 8.4b Adapted grid network at 6000 time steps (90 x 60) 94 Figure 8.6b Adapted grid network for converged solution for Mach 0.91 Figure 8.7a Comparison of computed solutions on fixed and adapted grid networks with TVD scheme 26 In the original iterative procedure, the most recently updated values for each grid point location were used in updating the remaining points and thus the procedure resembles a Gauss-Seidel iteration. This approach, however, is not suitable for vector processing and thus cannot benefit from the increased processing speed. Therefore, another approach is used in which all the grid points along a coordinate line are updated in a Jacobi iterative fashion. The new grid point values along that coordinate are then used in updating the points along the next coordinate line. Thus the technique can be viewed as a half Jacobi, half Gauss-Seidel iteration. This technique, which is suitable for vector processing, reduced the computer processing time for each iteration by 40 percent. Furthermore, only a few more iterations were required for the same convergence criterion in comparison to the Gauss-Seidel iteration, thus the effective speed up in processing time is approximately 40 percent. Ill.2 Modified Solution Procedure As pointed out earlier, the grids used for viscous transonic flow problems contain highly refined grid spacing in the direction normal to the surface (rÂ¡ direction) in order to resolve the viscous sublayer. This spacing results in grid cells with aspect ratios, up to 105 in the present applications, near solid boundaries. Iterative solution procedures for elliptic equations become inefficient for such large aspect ratios since the motion in one direction will be limited by the small spacing in the other direction. The effect of the aspect ratio on an iterative process can be ascertained by considering the simple grid cell shown in Figure 3.1 which for convenience is rectangular and aligned with the Cartesian 72 Figure 6.3c Grid network with internal boundaries Figure 6.4 Grid points near internal boundary 59 Figure 5.10 Initial fixed grid network for projectile base flow problem (70 x 60) 106 12. Ames, W.F., "Numerical Methods for Partial Differential Equations.11 Academic Press, Inc., New York, 1977. 13. Tu, C.G., "Development of An Adaptive Grid Generation Code for Projectile Aerodynamic Computation," Ph.D. Dissertation, University of Florida, 1985. 14. Pierson, B.L. and Kutler, P., "Optimal Nodal Point Distribution for Improved Accuracy in Computational Fluid Dynamics," AIAA Journal, Vol. 18, NO. 1, pp. 49-54, 1980. 15. Thompson, J. F., Warsi, Z. U. A. and Mast in, C. W. Numerical Grid Generation: Foundations and Applications. North-Hoiland, New York, pp. 215-221, 1985. 16. Nietubicz, C.J., Heavey, K.R. and Steger, J. L., "Grid Generation Techniques for Projectile Configurations," Proc. 1982 Army Numerical Analysis and Computers Conference, ARO Rept. 82-3, February 1982. 17. Warming, R.F. and Beam, R., "On the Construction and Application of Implicit Factored Schemes for Conservation Laws," SIAM-AMS Proceedings, Vol. 2, pp. 85-99, 1978. 18. Pulliam, T.H., and Steger, J.L., "Implicit Finite-Difference Simulations of Three-Dimensional Compressible Flow," AIAA Journal, Vol. 18, NO. 2, pp. 167-169, 1980. 19. Nietubicz, C.J., Pulliam, T.H. and Steger, J.L., "Numerical Solutions of the Azimuthal Invariant Thin Layer Navier-Stokes Equations," AIAA Paper 79-0010, Presented at AIAA 17th Aerospace Sciences Meeting, New Orleans, Louisiana. 20. Anderson, D.A., Tannehill, J.C. and Pletcher, R.H., Computational Fluid Mechanics and Heat Transfer. McGraw-Hill, New York, p. 546, 1984. 21. Baldwin, B.S., and Lomax, H., "Thin Layer Approximation and Algebraic Model for Separated Turbulent Flows," AIAA Paper 78-257, Presented at AIAA 16th Aerospace Sciences Meeting, Huntsville, Alabama, 1978. 22. Yee, H.C., "Linearized Form of Implicit TVD Schemes for the Multidimensional Euler and Navier-Stokes Equations," Comp. & Math, with Appls., Vol. 12a, Nos. 4/5, pp. 413-412, 1986. 23. Shiau, N. and Hsu, C.C., "A Diagonalized TVD Scheme for Turbulent Transonic Projectile Aerodynamic Calulations," submitted to 8th Int. Conf. on Computing Methods in Applied Sciences and Engineering, Versailles, France. 24. Sahu, J., Nietubicz, C.J. and Steger, J.L., "Numerical Computation of Base Flow for a Projectile at Transonic Speeds," ARERL-TR-02459, U.S. Army Ballistic Research Laboratory, June 1983. 98 the base region is shown in Figure 8.9b. There is generally good agreement between the two calculated solutions with the experimental data over the projectile surface with a slight discrepency near the shock on the cylinder. At the base comer and along the base the two calculated solutions are different though their trends are similar. Without any experimental data for comparison it is difficult to judge the accuracy. Furthermore, the thin layer approximation and the current turbulence model may not be valid near the base since a separated wake region exists. The adaptive grid network corresponding to the converged solution is shown in Figure 8.10a and an enlarged view of the grid network near the base comer is shown in Figure 8.10b. In comparing this grid network to that of Figure 8.5c for the projectile with sting problem it is evident that the r? coordinate lines emanating from the front portion of the projectile do bend upstream, but not so drastically, and do not prevent the convergence of the solution using the TVD scheme. There is the same general adaptation to the expansion waves and shocks, however the adaptation to the shock over the boattail appears to bend downstream, a result due to the smoothing of the control functions in the computational plane. Adaptation of the grid network near the base comer region is also evident. Note that the grid lines remained nearly orthogonal in the base region even though the normal direction changes direction rapidly near the sharp comer. In future work, the orthogonlity near the comer may be improved by decreasing the spacing of the Â£ coordinate lines in the reduced grid network. Also, the grid spacing between the rÂ¡ coordinate lines is fairly uniform and smooth and the coordinate lines did not 64 the adaptive grid generation equations for the base flow grid network. Recall that thel control function Q is defined so that the original grid point spacing in the normal direction would be retained. However, the point spacing resulting as a solution to the adaptive grid equations differs along each rÂ¡ coordinate line as shown in Figure 6.2a. The attraction of points is most prominent near the sharp comer at the projectile base where the effect is large due to the high curvature at the comer. A simple and effective technique to eliminate the effect of curvature can be developed by considering the solution of Laplace's equations between two concentric circles and C2 with radii R]_ and R2 respectively. V2Â£ = 0 (6.2.1) VrÂ¡ = 0 Solving Laplace's equation is equivalent to solving the adaptive grid equations with constant control functions and neglecting orthogonality (i.e. A = 0.0). In using constant control functions, we expect the coordinates to be equally spaced in each direction. It is advantageous to write the equations in polar coordinates (r,0), and it is sufficient here to seek a solution in which q is a function of 0 only and rÂ¡ is a function of r only. The Â£ coordinate lines then align themselves with the 0 coordinate direction and the rj coordinate lines became aligned with the radial coordinate r. The problem can then be formulated as 33 which, after using forward difference expressions to indicate the spacing between two adjacent points can be written ASj = AsjCj (3.2.9) C j + Qj+iQj-i Qj+i~Qj-i 2 (3.2.10) Now, as shown in Figure 3.5, if the ith point is removed it will be necessary to adjust the control function such that the new relationship is satisfied Asj + = ASj_! (3.2.11) where Cj_x is the modified control function and ASj_x is the increased grid point spacing ASj A Sj 4- ASj (3.2.12) The correct form of Cj_x can be obtained in the following derivation "k k Asj + 1 = ASj_x Cj _x ASj Cj + x (ASj_x+Asj) Cj_i Asj 1 Cj Cj+i = Asj _x (1+Cj) Cj_x C1C1 + 1 * = C* 1+C-i J 1 (3.2.13) 14 measures the orthogonality of the Â£ and rÂ¡ coordinates. As will be shown later, grid smoothness can be controlled through the definitions of the control functions. These three integrals can be combined into one total functional lip + hullL + a (VÂ£ .V7) 2 \ dxdy (2.2.4) p Q J where A is a parameter that weighs the relative importance of adaptation and orthogonality. Before applying the Euler-Lagrange equations for f and rÂ¡, it is beneficial to scale each term in the equations so that they remain of the same order of magnitude. An ordering analysis shows that the first term in eq. (2.2.4) is of order (C2/P), the second is of order (C2/Q) and the orthogonality term is of order (C4/!,2) where C is a computational length scale and L is a physical length scale. It is common practice to use highly stretched grids in transonic flow problems to provide extremely small grid spacing in the viscous sublayer region. The spacing along the normal coordinate can increase by five orders of magnitude, consequently the terms in eq. (2.2.4) may vary in a similar way. It is therefore beneficial to scale the terms locally rather than globally since a global scale cannot reflect the size of the term throughout the domain. The local scales Ap, Aq, and A0 used here are Ap=pL Aq=QL ao=jL (2.2.5) where Pp and Qp are the local values of the control functions and Jl is the local Jacobian and is of order (L2/C2). Although they vary throughout the domain they are considered constants during the application of the 83 Q do6) Figure 8.1 Comparison of control functions Figure 8.2 Lack of convergence using a self-adaptive computational technique with Beam and Warming scheme 73 be identical at the two points adjacent to the internal boundary. Referring to Figure 6.4, if the h coordinate line designated as an internal boundary contains the i^*1 point along a Â£ coordinate line, we then require ^ . se . i-1 . se . i+1 (6.3.1) To implement this requirement into the solution procedure, consider the one dimensional equation for adaptation along a Â£ coordinate, eq. (2.3.8) which is written here as P fii (6.3.2) The left hand side of eq. (6.3.2) appears in the source terms of the two dimensional equations for adaptation, eqs. (2.3.7). Substitution of eq. (6.3.2)into the source terms for a grid point just ahead of the internal boundary results in the modified source terms at the point i-1 (si)i-l,j ai-l,jJi-l,j The modified source terms for the point just behind the internal boundary are obtained by replacing i-1 with i+1 in eq. (6.3.3). In order to apply the iterative solution technique a value for the terms (s^/s^) in eqs. (6.3.3)for each point adjacent to the internal boundary will be needed. Their final values are not known a priori since they will depend on the final position of the points in each section, therefore they are set to (6.3.3) i-l,i 90 coordinates any change in the grid network used will change the dissipation terms and consequently alter slightly the transformed governing equations. Thus the continuous adaptation of the grid network continuously alters the dissipation terms and may be the reason the solution does not converge. VIII.3 Results Using the Self-Adaptive Computational Technique: The TVD Scheme Since the solutions obtained using the Beam and Warming scheme did not converge, another scheme, the TVD scheme was also used. Recall that the primary difference between the two schemes is the dissipation terms added to the governing equations. The initial grid network used is that of Figure 5.4 with IM=70 and JM=60. The time step used is 0.1 and the grid was adapted every 200 time steps. For Mach 0.96, the solution converged easily within 2000 time steps. Only one iteration was needed for the grid network to converge at 1800 and 2000 time steps which also indicates a converged solution. Figure 8.5a shows the computed pressure coefficient using the adapted grid generation procedure compared to the solution obtained on the fixed grid network of Figure 5.4. The solutions agree well with experimental data over most of the projectile surface, except for the predicted position of the shock over the cylinder. The position predicted using the adaptive grid generation technique is centered on the cylinder whereas the position predicted using the fixed grid network is slightly upstream of the center. Figure 8.5b shows a shadowgraph for corresponding flow conditions in which it can be seen that the experimentally predicted shock location is centered over the cylinder portion of the projectile. Thus, the use of the adaptive grid generation technique improved the solution accuracy. Figure 8.5c shows 99 overlap at the comer, a result attributed to the elliptic nature of the adaptive grid generation equations. The adaptive grid generation procedure, thus, can successfully provide a good adapted grid for this complex geometry. 107 25. Kayser, L.D. and Whiton, F., "Surface Pressure Measurements on a Boattailed Projectile Shape at Transonic Speeds," ARBRL-MR-03161, U.S. Army Ballistic Research laboratory, March, 1982. 20 The one dimensional equation for the boundary adaptation, eq. (2.2.9), is also transformed in a similar manner. The second derivative term can be shown, using the chain rule, to transform as sÂ£ Â£ ss = T (2.3.8) sr Substituting eq. (2.3.7) into eq. (2.2.9) yields the transformed one dimensional equation for adaptation along Â£ boundaries + sÂ£ PÂ£/P =0 (2.3.9) and for adaptation along boundaries coincident with rÂ¡ coordinates s^ + s^Q^/Q = 0 (2.3.10) These equations can be solved directly for the grid point spacing along a coordinate. For instance, integrating eq (2.3.10) once yields the grid point spacing s r, = c/q (2.3.11) where c is a constant of integration and can be shown to be C = (2.3.12) f da. J Q where Sp> is the total arclength along the coordinate in the physical domain. However, since the solution in the interior will require an iterative algorithm, the differential form, eqs. (2.3.9) and (2.3.10) are 2 be concentrated to provide high resolution in certain areas but will always form an equally spaced rectilinear coordinate system in the transformed plane. Thus, once the governing equations are transformed onto the curvilinear coordinates the finite difference algorithms used to solve these equations may be developed on the equally spaced rectilinear coordinate system which is inherent in the transformed plane. Certain characteristics of the general curvilinear coordinate system have been shown to affect the accuracy of the finite difference approximations in the transformed plane [3]. These include orthogonality of intersecting coordinate lines and the smoothness of the grid spacing. It is also known that a refined mesh is needed in regions of large physical gradients to obtain accurate approximations. In fact, poor resolution in these high gradient regions can be detrimental to solution accuracy and to the convergence of finite difference algorithms for fluid flow problems [4]. Thus, the development of a good adapted grid network is important for the successful application of finite difference methods to solve complex fluid dynamic problems. The self-adaptive grid generation technique has became an important area in computational fluid dynamics since it has been shown to provide good grid networks for the complex flow fields occurring in transonic and supersonic flows [5,6,7]. The use of boundary fitted curvilinear coordinate systems with transformed governing equations leads naturally to the concept of solution adaptive grids. As constraints of computer storage and CRJ time prohibit uniform mesh refinement throughout the entire domain, the coordinate spacing in the physical domain is varied to increase resolution in only those regions in which the solution is changing rapidly. For simple flow problems, when the position of frOT %Aux uA$x = Â£ zp i* + ^ ^ x = Mi ^ 5aJx + i^i! |x + ueiS )xUxz ii| Sx = >tXq Ze 2e er/((A)xxT^x + (x)xxT A-) = xxk Er/((A)^xq5x + (x)^xt A-) = ^XU Er/((A)^^qx + (x)^^q A-) = xxu er/( (A)xxq^x + (x)xxq **A-) = Er/( (AJ^q^x + (x)^xq UA~) = ^x5 cr/((AJ^^I^X + (x)AAq ^A-) = xx AXrj XIQNaddV 45 as the transonic flow considered here, large velocity gradients are expected to occur normal to the surface and most of the grid points are used to resolve these gradients. As larger spacing is then used along the surface and the viscous derivatives in that direction are expected to be small in comparison, they are dropped from the equations. The conservative form of the transformed thin layer governing equations for axisymmetric flow, written in the compact vector form are atq + dtk + dÂ£ + b = (s.i.i) where p P U A : pv , E=J* p w . C , P U pV puU+?xp A X pUV+r?xp pvU+^yp , F=J* pW+flyP pwU+ezp pWV+r7Zp (e+p)U j . (e+p)V 6=J* 0 0 -pw2y p/y pW(y^U+y^V) 0 (5.1.2) A s = J* 0 uRu^ + (u/3)x uRv,, + (u/3)T*y }, uRw^ + (u/3)Tr)z A(%m (u2+v2+w2 ) ^uPf1 (7-I) _1 (a2) ^ + (u/3) (r?xU+r?2yv+r?2xw)T and R = ^x2 + Vy2 + 1z2 T + nyVr, + *?xw>7 A SELF-ADAPTIVE COMPUTATIONAL METHOD FOR TRANSONIC TURBULENT PROJECTILE AERODYNAMICS BY CHRISTOPHER WILLIAM REED A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FIORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FIORIDA 1987 40 where 7p and 7q are parameters and f and g are some derivatives of the dependent variables scaled to range between zero and one. In the following developments, the terms Wf7 and h will be used to indicate either P,7p and f or Q,7q and g. By evaluating numerically the solution to the one dimensional equation for adaptation, the following expression equating the control function to the grid point spacing between points (i) and (i+1) along a coordinate line is obtained, ASi Sijr 1 l+7hj S j 1 l+7hj (4.1.2) After evaluating this expression for the maximum and minimum spacing corresponding to h = 1 and h = 0 respectively, and forming their ratio (Asi) max = 1 + 7 (4.1.3) (Asi)mim it can be seen that the specification of 7 dictates the ratio of maximum and minimum spacing along a coordinate line. The parameter 7 also influences the smoothness of the resulting grid point distribution. If it is chosen as zero, equal spacing as indicated by eq. (4.1.2) would result; as 7 is increased, the spacing must change more rapidly to obtain the varied spacing. Following the work of Nakahashi and Deiwert [7], another parameter a has been introduced to obtain more control over the spacing. The general control function now becomes W 1 + 7ha (4.1.4) 37 are approximately equally spaced and the spacing along the r? coordinate lines in the initial reduced grid is retained which indicates that the modified control function is correct. Upon inserting the removed points along the rÂ¡ coordinates using eq. (2.3.12), the final grid network, Figure 3.8, is obtained, which is the same grid that would have been obtained in Figure 3.3 after a large number of iterations. Ill.3 The Adaptive Grid Generation Code A computer code base on the iterative solution procedure described in the preceding section has been developed and is referred to here as the Adaptive Grid code. As the solution process is iterative, an initial guess for the grid point locations is required and the flow variables for which the grid network is to be adapted cure needed for calculation of the control functions. These are considered as input to the Adaptive Grid code. The initial guess for the grid network can be a grid network obtained by any conventional grid generation technique or a previously adapted grid network. Values for the flow variables at the grid points are obtained using the flow solvers which are described in Chapter V. The sequence of calculations in the Adaptive Grid code proceeds as follows: 1. Input an initial grid network and the flow variables at each grid point. 2. Calculate the control functions P and Q. 3. Remove points to form the reduced grid network and calculate the modified control function Q*. 4. Solve the grid generation equations for the reduced grid network using P and Q*. 5. Reinsert the removed points along the coordinate lines using the control function Q which then results in the adapted grid network. 23 T1 = (xi+l,j+l + xi-l,j-l xi+l,j-l ~ xi-l,j+l)/4 (3.1.Id) T2 = (Yi+1,j+1 + Yi-1,j-1 Yi+1, j-1 Yi-1,j+l)/4 where the source terms and the coefficients all contain first derivatives which vhen approximated with central differences do not contain the point xi,j' Yi,j explicitly. This set of equations is solved using a Newton- Raphson iterative scheme. The initial guess for the grid point locations will not satisfy the finite difference representation of the equations and therefore a non-zero residue will exist (Rl)l,j = eq (3.1.1a) (3.1.2) (R2)i,j = eq (3.1.1b) where i. indicates the iteration number. An estimate of the residue at the next iteration can be obtained by expanding the residue in a Taylor series about the current residue and retaining only the linear terms, 3(Rl)i,j a^i,j l a 0 4-1 o a cd- \ . (3.1.3) Â¡L d (r2)i,j aYi,j where ax and Ay are the changes in the grid point location 5 Saltzman and Brackbill [5] have developed an adaptive grid generation scheme based on a variational approach. A functional is defined to measure each grid characteristic and the minimization of these functionals results in a set of partial differential equations that govern the grid point spacing. With the proper choice of parameters the equations are elliptic, which guarantees a one-to-one mapping and results in a smooth grid. They solved an inviscid supersonic flow problem by adapting the grid cell size to a function of the pressure gradient to capture shocks. However, it is not possible in this approach to adapt the grid spacing independently in each coordinate direction, a feature that is prohibitive if the spacing in different coordinate directions vary by a large amount. The adaptive grid generation technique presented here is similar in approach to that of Brackbill and Saltzman but is developed for applications to viscous transonic flow problems. The solution of these problems usually contains shock structures aligned with one coordinate direction and boundary shear layers parallel to a streamwise coordinate. Also, the grid point spacing can vary by orders of magnitude along different coordinates and it is thus necessary to adapt the grid independently in each coordinate direction. Therefore, a functional is defined for independent adaptation in each coordinate direction such that the resulting equations are elliptic. The use of elliptic equations in grid generation offer some important advantages that motivate their use. They constitute a boundary value problem and can thus be used in any domain. They guarantee a one to one mapping and also, grid smoothness is inherent in the equations. In addition to the equations governing the grid, there are two other topics that require attention when developing 68 P =( Y V*r,r, a 3/2 (6.2.9) By multiplying this result with the inverse of the spacing along a Â£ coordinate (l//y), the following expression is obtained Ki = Yr)^-r)r] ~ ^-r)Yr]r] Jy 3/2 (6.2.10) which is the term that is subtracted from the adaptive grid equation for Â£ in order to cancel the effects of curved rÂ¡ coordinate lines. An expression for the influence of curved Â£ coordinates on the spacing along rÂ¡ coordinate lines can be obtained in a similar manner and is K2 = y^x^ xcyee Ja 7 3/2 (6.2.11) Thompson et al. [15] have derived similar expression, but do not implement them in the same way. It should be noted that in the computational plane, the terms and K2 contain second derivatives of the variables x and y, and when eqs. (2.2.7) are transformed onto the computational plane, the second order central difference approximations to these terms depend explicitly on the point (xifj, yi,j) However, in implementing the Newton-Raphson iterative solution procedure, the terms Ki and K2 are considered part of the source terms. Thus, the right hand side of eqs. (2.3.7) becomes CHAPTER III NUMERICAL SOLUTION OF THE ADAPTIVE GRID GENERATION EQUATIONS III.l The Newton-Raphson Iterative Method Equations (2.3.3), (2.3.9) and (2.3.10) constitute a set of nonlinear elliptic partial differential equations. In order to solve these equations numerically they are approximated with second order central difference expressions which results in a set of coupled algebraic equations for each grid point. The finite difference approximation for the equations at grid point, (Xj^j, Yifj) become al(xi+l,j 2xi,j + xi-l,j) + a3(xi,j+1 "2xi,j + xi,j-l) + bl(Yi+l,j -2yij +Yi-l,j) + b3(Yi,j+1 -2Yi,j + 2Yi,j-l)-Fi=0 (3.1.1a) cl(xi+l,j "2xi,j + xi-l,j) + c3(xi#j+1 -2xi/j +Xifj_i) + dl(Yi+i,j -2Yi,j +Yi-l,j) + d3 (Yi, j+1 -2Yi,j +Yi, j-l) -F2=:0 (3.1.1b) P|Ja Fi = a2Ti b2T2 (3.1.1c) Qr;J7 F2 = a2Tx d2T2 22 55 MACH 0.96 cPt i Figure 5.6 Converging solution for the Beam and Warming scheme Figure 5.7 Converged solution for the Beam and Warming scheme ACKNOWLEDGEMENTS The support of Dr. Chen-Chi Hsu has been most encouraging during the course of this work. His guidance and willingness to provide it are appreciated. In giving me the freedom to experiment, discover, fail and succeed I have learned as much about research as I have about the topic. For this I am grateful. My appreciation is extended as well to the committee members, Dr. Leadon, Dr. Kurzweg, Dr. Mikolaitis and Dr. Wilson. This work was partially supported by the APOSR. Computer time was provided by NASA Ames, NSF and the Pittsburgh Supercomputer Center. iii 101 Figure 8.10 Adapted grid network for the converged solution for Mach 0.96 96 The results for a third case, Mach 1.2, are shewn in Figures 8.7a and 8.7b. There is no noticeable difference between the computed results on the fixed grid network and the results obtained using the adaptive grid generation procedure except near the projectile nose where the grid configuration differs the most. In the adapted grid network, clustering again is evident near the junctures, however no shocks are present for this case, a result also evident in the grid network. The equal accuracy of the solutions obtained with the adapted grid generation procedure and on the fixed grid network is a result of the similarity of the grid networks. In comparing the adapted grid networks of Figures 8.6b and 8.7b to the fixed grid network of Figure 5.4 it is evident that they have basically the same features in that the grid point clustering along the Â£ coordinates is concentrated near the two junctures. Thus the original fixed grid network is well adapted to the problem since the shocks for case of Mach 0.91 are so close to the expansion waves and no shocks occur in the case of Mach 1.2. However, it should be noted that the adaptive grid generation procedure successfully produced a good adapted grid network without any knowledge of the position or orientation of the important flew phenomenon. The choice and construction of the proper control functions of course are important in obtaining these results. In the context of these results it can be concluded that the adaptive grid generation procedure can reliably produce a good adapted grid network provided good choices are made for the control functions. 82 where the derivatives are with respect to the arclength s along each Â£ coordinate line. This function is largest at the expansion waves where the pressure gradient changed sign. At the base and peak of a shock eq. (8.1.1) also increased but it is not as large as the value near the expansion waves. Thus the smallest spacing will occur near the expansions and yet the grid points will also cluster, to a lesser extent, in the vicinity of shocks. In the normal direction, the important area requiring increased resolution is the boundary layer and the logical choice for the function g in the control function Q is the velocity gradient. However, it was found that the control function based on the velocity gradient resulted in too many points being placed in the boundary layer and led thus to a rapid stretching of the grid point spacing along the ry coordinate lines. Another choice for the function g, the exponential of the velocity gradient g = exp dn d s (8.1.2) where u is the velocity magnitude and s measures along the rÂ¡ coordinate line, was found to yield better results. Figure 8.1 shows a plot of the control function Q along a typical Sv coordinate line. The control function resulting from the exponential of the velocity gradient is similar to the control function corresponding to the grid point distribution in the fixed grid networks which is known to give good results provided enough points are used. The control function obtained using the velocity gradient however is large for many grid points and, 63 Figure 6.1a Grid network using X=0.0 Figure 6.1b Grid network using \ =0.5 57 shocks. However, after approximately 5000 time steps the solution does converge and agrees well with the experimental data as indicated by the comparison of the experimental and computed pressure coefficient along the projectile surface. Figure 5.8 shows the results obtained using the TVD code. The solution also agrees well with the experimental data; moreover it converged within 1600 time steps. This difference in the convergence ofthe two schemes is shown in Figure 5.9 which plots the least square residue R, versus the number of time steps N. Here the slower, oscillatory convergence of the Beam and Warming scheme is evident while the TVD scheme appears to converge monotonically. The projectile base flow problem is also solved on a fixed grid using the TVD code. The grid network for this calculation was obtained using a modified version of the technique of Steger et al. [16] in which the boundary condition on the downstream bounding rÂ¡ coordinate line has been changed so that the Â£ coordinate lines intersect with the x axis at right angles. The grid network obtained with IM=70 and JM=60 is shown in Figure 5.10. As can be seen, the Â£ coordinate lines now bend around the sharp comer, and end at the downstream axis of symmetry forming an O-type grid network. The computed pressure coefficient obtained for this problem using the TVD code is shown in Figure 5.11. Again, the solution converged in approximately 1600 time steps. It compares well, as before, over the secant-ogive and cylinder portions but then differs from the experimental results at the end of the boattail. No experimental data is available for the base region. There are a number of reasons that may cause the difference between the computed and experimental results near the comer. One reason may be the presence of a narrow sting used to support the experimental model that is not present in the computational 36 Figure 3.7 Reduced grid network after 20 iterations Figure 3.8 Final grid network after reinserting grid points 80 the initial grid. For the purpose of adapting the grid network, the source terras corrected to eliminate curvature effects, eqs. (6.2.12) were used in the adaptive grid generation equations. In solving the projectile problem with sting two internal r? coordinates were designated as internal boundaries. The first is over the secant ogive as shown in Figure 6.5, the second separates the grid network over the projectile from that over the sting and is located at approximately x=7.0. In solving the projectile base flow problem, no internal boundaries were used. To Shelley and Drew 21 retained and solved iteratively so that the boundary points remain consistent with the interior points as the grid evolves. CHAPTER IX CONCLUDING REMARKS The self-adaptive computational technique developed here has been shown to provide good adaptive grid networks for the transonic projectile flow problems considered. The equations governing the grid adaptation were developed, based on a variational approach, to provide adaptation independently in each coordinate direction and also enhance both grid orthogonality and smoothness. The resulting set of elliptic partial differential equations provide explicit control of orthogonality and adaptation while smoothness is inherent in the ellipticity of the governing equations. Smoothness can also be controlled in the definition of the control functions. A series of modifications were made to improve both the adaptive grid generation equations and their solution procedure. local scaling of the equations was found useful for highly stretched meshes, and an effective technique to eliminate the effects of curved boundaries, especially useful in the projectile base flow problem, was developed. Also, the numerical solution procedure for the adaptive grid generation equations was improved to increase the poor convergence rate due to the high aspect ratio grid cells present in the grid networks used for transonic viscous flow problems. Finally, the use of internal grid boundaries, which were used in the projectile flow problem with sting, helped to provide a better grid network by restricting the movement of the grid points. 102 75 Figure 6.5 Smoothed grid spacing near the internal boundary 30 coordinate will be retained in the solution to the adaptive grid generation equations. Using the Newton-Raphson iterative solution procedure, the solution is advanced 20 iterations which results in the grid shown in Figure 3.3. As can be seen, the point spacing along the Â£ coordinates is approximately equal for the coordinates sufficiently above the surface. However, in the region of the high aspect ratio cells, the movement of the points in the Â£ direction is retarded by the extremely small spacing in the x] direction. A modification has been made to eq. (3.1.11) which updates the points along the bounding Â£ coordinate in formulating this example. The movement of points lasing the one dimensional equation is independent of the interior points and thus is not restricted by the high aspect ratio cells. This creates an inconsistency at the boundary since the points adjacent to the boundary are constrained by the high aspect ratio cells. Equation (3.1.11) therefore has been modified as i a i ASi = (Rsi)i ( ) (3.2.6) 1+AR. l where AR* is the aspect ratio of the cell adjacent to the 1th point along the boundary Â£ coordinate. This change only effects the rate of convergence of the iterative solution procedure for the boundary points so that their movement will remain consistent with that of the interior points adjacent to the boundary during the convergence process. In order to alleviate the deficiency of the above solution procedure, a technique has been developed which is based on solving for a reduced grid in which many of the points in the r? coordinate direction have been 51 Another problem considered is the projectile base flew problem. In Figure 5.3, which shows the basic O-type grid configuration, the Â£ coordinate lines bend around the sharp comer at the base and end on a downstream axis of symmetry which is now an rÂ¡ coordinate. The existing codes are used for this problem by simply modifying same of the boundary conditions. The boundary conditions along the upstream axis of symmetry, line AD, and the projectile surface, line CD, are the same as those used in the projectile problem with sting. For the downstream axis of symmetry, line BC, the same conditions along the upstream axis are used; however, a backward finite difference formula is required. The freestream conditions are applied along the outer boundary along the line AA1 and the downstream boundary conditions for the projectile with sting problem are applied along the portion of the outer boundary between A'B, except that the extrapolation is done along the rÂ¡ coordinate lines. V.3 Examples on Fixed Grid Networks In order to show the general characteristics of the Beam and Warming and TVD schemes, and for purposes of comparison with the solutions obtained with the adaptive grid generation technique, the two schemes were used to solve the projectile flow problem with sting on a fixed grid network. The fixed grid network was obtained using a method developed by Steger et al. [16] which is an effective technique for external flow problems. Two equations, one enforcing orthogonality and one specifying the grid cell volume, form a set of hyperbolic partial differential equations which are solved numerically to generate the grid network. In the marching solution algorithm, the initial grid points are specified along the projectile surface. Ihe solution is marched outward in the v 31 Figure 3.3 Grid network after 20 iterations f Figure 3.4 Reduced aspect ratio cell 71 the secant portion of the projectile. The solution did not converge and instead the oscillations continued and eventually effected the solution downstream. In order to prevent the upstream bending of the r? coordinates, one rÂ¡ coordinate line in the initial grid network is designated as an internal boundary. This rÂ¡ coordinate line remains fixed in space during the iterativesolution procedure, thus no other rÂ¡ coordinate lines will pass through it. However, the points defining the internal boundary coordinate line are allowed to move along the coordinate so that they remain consistent with the spacing of the points along the adjacent rÂ¡ coordinate lines. This is accomplished by updating the arclength distribution along the internal boundary after each iterative sweep in which the new distribution is defined as the average of the arclength distributions along the two adjacent rÂ¡ coordinate lines. Figure 6.3c shews the grid obtained when the 12th rÂ¡ coordinate line from the projectile nose is designated as an internal boundary and as can be seen it prevents the rÂ¡ coordinate lines from bending upstream. It has been chosen to lie in a region in which no major grid clustering is expected so that it will not effect the adaptation of the grid network. The length and number of points used along the Â£ coordinate lines in each section may result in an abrupt change in the point spacing along the Â£ coordinate lines passing through the boundary as can be seen in Figure 6.3c. Therefore the control function P, controlling the spacing along the Â£ coordinate lines is locally modified near the internal boundary so that the grid point spacing varied smoothly through the internal boundary. This can be accomplished by requiring the normalized rate of change of grid point spacing along a Â£ coordinate line (s^/s^) 3 important solution gradients is known, good adapted grids can be obtained with conventional techniques. However, in more complex problems the position of these important regions is not known a priori and then the development of a proper adaptive grid network is difficult. Solution adaptive grid generation adresses this problem by continuously updating the grid network as the solution evolves such that the important physical gradients are sufficiently resolved as they develop. It is thus an attempt to efficiently reduce the truncation error by only refining the mesh in those regions where the error may be large. Also, since the grid characteristics of smoothness and orthogonality also affect the truncation error, a good adaptive grid generation technique should include the enhancement of these characteristics as well. The general approach of most adaptive grid generation schemes is based on minimization techniques. A measure of each desired grid characteristic is defined, and the governing equations for the grid are obtained by minimizing the integral of these measures over the domain. In many instances adaptation is required in only one coordinate direction. Ewyer et al. [8], for instance, have developed an equidistribution law to control the grid spacing along one coordinate direction. The spacing was adapted to a flame front in a two-dimensional combustion problem and the improved resolution eliminated spacial oscillations in the flame front. Gnoffo [9] has developed another one dimensional adaptive scheme in which tension springs are assumed to connect the grid points along a coordinate. The spring constants became large when increased resolution was required, and the points clustered in those regions in an attempt to minimize the potential energy of the spring system. The application of this technique to supersonic flow over a Viking Aeroshell successfully 32 removed to decrease the maximum cell aspect ratio. Figure 3.4 shows a diagram of a typical grid for viscous flow problems containing large aspect ratio cells near a solid boundary. A reduced grid can be formed by removing all the Â£ coordinate lines denoted with the dashed lines, resulting in a grid with fewer points along the r? coordinate lines which are designated by the intersection of the solid lines. Currently, enough points are removed along the rÂ¡ coordinate lines such that the minimum spacing between Â£ coordinate lines is greater than 0.04. The typical spacing along the boundary is of order (0.1) which results in cell aspect ratios of order (1). Upon convergence of the solution for the reduced grid, all the points removed are reinserted along coordinate lines using the same one dimensional equation that controls the spacing along the boundaries. It is necessary in this process to adjust the control function Q at the remaining points to account for the spacing requirements of the points removed. A simple procedure for this purpose can be derived by using the equation for one dimensional adaptation, eq. (2.3.10). The procedure will also be adequate for the two dimensional adaptive grid generation equations. If we write explicitly the finite difference approximation for eq (2.3.10), sj+l-sj-l Qj+lQj-l sj+l-2sj + sj_! + ( )( )/Qj=0 (3.2.7) we can regroup the terms in the form Qj+lQj-l Qj + ( ) (sj+l-sj)=(Sj-Sj_1) ( ) (3.2.8) Qj+l"Qj-l Qj ( ) 49 where D and V are forward and backward difference operators, respectively. Since the scheme is not naturally dissipative these terms are added to the finite difference scheme to prevent nonlinear instabilities and high frequency spatial oscillations. The constant coefficients ej and eg used here are 4At and 2At respectively which are consistent with values used in other applications of this code [19]. In general, these coefficients should be chosen large enough to dampen the instabilities but kept small enough to prevent any degradation of solution accuracy. In Chapter VII, seme examples using the Thin Layer code with the Adaptive Grid code shew that there are problems in obtaining a converged solution for the thin layer equations when the grid is adapted. Therefore another numerical scheme for solving the thin layer equations is also used. This scheme is the TVD scheme developed by Yee [22] for the multidimensional Navier-Stokes equations. The theory and development of TVD schemes is beyond the scope of current discussions. However, for purposes here it is sufficient to view the TVD scheme as an improvement in the artificial damping terms in the Beam and Warming scheme. In fact, the TVD scheme developed by Yee for multidimensions can be implemented into the existing Thin Layer code based on the Beam and Wanning scheme by modifying only the added dissipation terms. The original Thin-Layer code has been modified to provide an option for the TVD scheme [23] and is referred to here as the TVD code. These 'sophisticated' dissipation terms switch the scheme from second order to first order accuracy at points of extrema and result in a more dissipative scheme in those regions. It should be pointed out that the primary purpose of the current research is to develop an adaptive grid generation technique. The TVD scheme is used 41 By prescribing the minimum spacing and 7, eq (4.1.2) can be written as (Asi)min ? (i+7h? ^ 1+7 = 0 (4.1.5) in which a is the unknown quantity and can be determined by a root finding technique. Currently, the Newton-Raphson finding method is employed. Thus, both the minimum and maximum grid point spacing along a coordinate can be specified by specifying the two quantities, (As^J^^ and 7. It should be noted that the prescribed values may not be realized in practice due to the coupling effects of the two dimensional adaptive grid generation equations. Also, the constraint of orthogonality may alter the results in some areas. IV.2 Smoothing The choices for the functions f and g will of course depend on the problem being solved. However, the functions will in general, contain derivatives of the solution which are approximated numerically and it is therefore useful to smooth the control functions to eliminate any spurious data. Also, since smoothing is most effective when the function changes most rapidly, it will help reduce any rapid variations in the control function and subsequently smooth the grid point spacing. The method of smoothing used here, which is similar to that used in reference [5] is 66 d2Â£ = O de2 d2r? 1 dr? + =0 dr2 r dr (6.2.2) Â£ = 6 on Ci and C2 rÂ¡ = R]_ on C-]_ 1 = R2 on C2 The boundary conditions corresponds to uniform spacing along the Â£ coordinate lines coincident with the bounding concentric circles. The spacing along each coordinate direction is the inverse of the first derivatives and rÂ¡r which by integrating eqs. (6.2.2) are found to be to ~ a *?r (6.2.3) where a and b are constants of integration. These results indicate that the spacing in the Â£ direction is uniform while the spacing in the rÂ¡ direction is smaller near the inner circle and larger near the outer circle C2 The term in the governing equation responsible for this result, designated here as K, is K 1 d^ r dr (6.2.4) which can be viewed as the curvature of the intersecting Â£ coordinate line, (1/r), multiplied by the inverse of the spacing along the rÂ¡ transonic projectile flow problem, the grid network was adapted to the pressure distribution in the streamwise direction and to the velocity distribution in the direction normal to the projectile surface. Results obtained for Mach 0.91, 0.96 and 1.2 indicate that the procedure can provide good adapted grid networks provided good choices are made for the control functions. vii 48 Layer code. The code employs Sutherland's law [20] to relate the viscosity and heat conductivity to the temperature and the effects of turbulence are modeled using the eddy viscosity model of Baldwin and Lomax [21]. The numerical algorithm is a time-marching scheme in which steady state solutions are obtained as the time asymptotic solution of the unsteady equations. Spatial derivatives are approximated with central differences and a forward difference approximation is used for the unsteady terms. The time marching scheme of the finite difference approximation to the governing equations are written in delta form as (I + h S^k ciDi) (I + hfi^ hR^-jfoj1 jD2 ) Aq = h (d^fe + 6n$ ~ Ri1 S fj 6) (5.2.1) e Ed3 where A and B are the Jacobian matrices a = a aÂ£ B = (5.2.2) aq and M is obtained as the time linearization of S. Here, 6 is the central difference operator. Details of this scheme can be found in Warming and Beam's [17] work, as well as in that of others. The terms D3, and D2 are implicit dissipation operators and D3 is an explicit dissipation term Di = J-1 VÂ£ J D2 = J_1 A J D3 = J"1 ((A^ Ve)2 + Ay V,)2) J Â§ (5.2.3) 103 In coupling the adaptive grid generation procedure with the flow solvers it was shown that the adaptive grid generation equations could provide good adapted grid networks for the transonic projectile flow problems. It was found that the control function based on the curvature of the pressure distribution was more effective than the pressure gradient when both shocks and expansion waves are present. Also, in using a control function based on the velocity gradient, it was shown that the adaptive grid generation procedure could provide the extremely refined grid point spacing necessary in the boundary layer region. The most important feature of the adaptive grid generation procedure developed here is the simple relationship between the control functions and the grid point spacing which is maintained in an elliptic set of partial differential equations. As more complex problems are attempted, the control functions for each coordinate direction can be carefully designed, based on various physical gradients, prescribed grid point distributions and combinations thereof, to obtain the desired adaptation. The elliptic nature of the adaptive grid generation equations will be helpful, especially for complex geometries, in maintaining a one to one mapping and preventing highly distorted grid networks from occuring. 60 Figure 5.11 Converged solution for projectile base flow problem using the TVD scheme 65 mw ui/i Figure 6.2a Grid network obtained without corrections for curvature Figure 6.2b Grid network obtained with corrections for curvature 100 0. 0. -0. -0. 2 0 a 2 4 / a - MACH 0.96 b x s X\ f .-..V. - adapted grid fixed grid b i 0.37 0.0 Figure 8.9 Comparison of computed solutions for the projectile base flow problem with TVD scheme 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 Fhilosophy. Chen-Chi Hsu, Chairman Professor of Engineering Sciences 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 (fegree of Doctor of Philosophy. Bernard tadori Professor of Engineering Sciences 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. Z//tn cl ?/ ')<'<* Ulrich H. Kurzweg Professor of Engineering Sciences 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 Fhilosophy. David W. Mikolaitis Assistant Professor of Engineering Sciences 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. C* ll^bz David C. Wilson Associate 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. August 1987 liaAcjf &. Dean, College of Engineering Dean, Graduate School 56 Figure 5.8 Converged solution for the TVD scheme Figure 5.9 Least squared residue 12 in the computational domain T ~ "gX| fxxx2 xÂ£Â£ fxx (2.1.6) The first term is the usual truncation error term containing the grid point spacing x^ and the third derivative of the function. The second term is due to the varied spacing along the Â£ coordinate lines in the physical domain. It is thus possible to increase the truncation error if the coordinates spacing is not smooth. They have also shown that truncation error for a first derivative varies inversely as the sine of the angle between the curvilinear coordinates and therefore a lack of orthogonality can also increase the truncation error. Orthogonality is also important near solid boundaries for turbulent flow problems because turbulence models often are based on information in the normal direction. In the context of a coordinate transformation, the adaptive gridding can be viewed as an attempt to reduce the truncation error by adjusting the transformation metrics. This is accomplished by judiciously defining the curvilinear coordinates such that they vary smoothly, are nearly orthogonal, and concentrate in regions where the solution is changing rapidly. II.2 A Variational Approach The choice of a method to generate the curvilinear coordinates is arbitrary in that it has no dependence on the partial differential equations that are to be solved. Of the many methods available, however, the variational approach is attractive for purposes of adaptive grid generation because it can provide a simple, explicit means of incorporating desired properties into a grid generation scheme and can 92 Figure 8.5a Comparison of computed solutions on fixed and adaptive grid networks with TVD scheme F Figure 8.5b SOCBT shadowgraph for Mach 0.96 (reprinted from reference 24) 86 Figure 8.3a Computed solution between 4000 and 5000 time steps with Beam and Warming scheme Figure 8.3b Computed solution between 5000 and 6000 time steps with Beam and Warming scheme CHAPTER I INTRODUCTION The study of boundary fitted grid systems has became an important part of computational fluid dynamics and, in fact, their development has been cited as a major pacing item in the progress of computational fluid dynamics [1]. The use of such systems makes the application of finite difference techniques versatile and effective in solving complex fluid dynamic problems in arbitrary domains. Boundary fitted grid generation can be viewed as the development of a general curvilinear coordinate system in which the coordinates coincide with the boundaries on the edges of the physical domain. This feature of the boundary fitted grid systems facilitates the application of finite difference approximations at the boundaries and eliminates the need for interpolation between grid points near the boundary. This is particularly advantageous when the boundaries are highly curved or have slope discontinuities since the use of interpolation near these boundaries may have a significant effect on the solution [2]. In the domain interior, the intersection of each family of coordinate lines denotes each grid point and therefore defines the computational mesh. The distribution of the coordinate lines in the physical domain may 1 34 Once Cj_i is calculated losing eq. (3.2.13), the modified value of Qj-i * can be obtained by solving eq. (3.2.10) for Qj-i 1+C1-1 Q*-l = ( ) (Qj + i-Qj-2) (3.2.14) 1-Cj_i This procedure is applied for each point along an coordinate that is removed. It should be noted that this procedure does not affect the final grid, the points remaining in the reduced grid will obtain the same location as if all the points were used in the iterative solution procedure. The removed points are, once the solution for the reduced grid is obtained, inserted back along the r? coordinates lines. This is done by approximating the r? coordinate with straight line segments between the points of the reduced grid and then reinserting the removed points along each segment with spacing dictated by eqs. (2.3.10) using the original control function Q. If necessary, the use of a higher order approximation to the r? coordinate lines, such as a cubic spline, can be used to obtain a better representation of the rÂ¡ coordinate lines. In using this procedure with the two dimensional grid generation equations, the solution may vary somewhat since the equations are coupled. Also, the effects of orthogonality may alter the results in some areas. In order to show the advantage of this procedure, the initial grid network of the previous example is advanced 20 iterations using the reduced grid technique. Figure 3.6 shows the initial reduced grid in which the points responsible for the high aspect ratio cells cure removed. After modifying the control function Q, the solution to the adaptive grid generation equations is advanced 20 iterations resulting in the grid shown in Figure 3.7. As can be seen, the points along all Â£ coordinates 85 attempts to use the adaptive grid technique, the grid network was adapted every 200 time steps. A series of calculated pressure coefficients obtained with this approach are shown in Figure 8.2. As can be seen, the pressure coefficient oscillates widely over the projectile surface and shows no indication of converging. In a series of attempts to obtain better results, the time step was decreased to 0.05, the dissipation terms were doubled and the number of time steps between adaptation were both increased to 500 and reduced to 50. However, similar results were obtained in each case. In the next case tried, the number of points in the streamwise direction was increased to 90 points, thus IM=90, and the grid network was adapted every time step. The results for this case are much better, however a converged solution was not obtained. Figure 8.3a shows the computed pressure coefficient at the times indicated, and as can be seen, the solution appears to have converged everywhere except in the shock boundary layer interaction regions. The shocks propagate upstream and downstream in a cyclic fashion. Figure 8.3b shows the same phenomenon occurring at a later time. In fact, the solution was continued for 8000 time steps and the propagation of the two shocks showed no indication of diminishing. Figure 8.4a shows the reduced grid network which resulted at 5000 time steps. As is evident, the grid point spacing in the rÂ¡ coordinate direction is clustered near the projectile surface with a distribution similar to that of the hyperbolic grid in Figure 5.4. The spacing varies somewhat in the shock and expansion wave regions, most likely as a result of enforcing orthogonality. The clustering in the streamwise direction CHAPTER VI PRELIMINARY TESTS AND MODIFICATIONS VI.1 The General Procedure Prior to coupling the Adaptive Grid code with the flow solvers to solve the projectile flow problem, it was found beneficial to modify the adaptive grid generation method. The following sections include an example showing the control over orthogonality and then describes two modifications made to the adaptive grid scheme to improve the general technique. It should be noted that each of the modifications is made to improve the method for the present flow problem and geometry and may not be necessary or useful in other problems. For each example the initial grid network used was that of Figure 5.5, or that of 5.10 if the base flow configuration is used in the example. In the adaptive grid code, the control function P was set equal to a constant so that the clustered points along each | coordinate line in the initial grid would tend to spread and equal spacing would result. The control function Q was defined using eqs. (4.3.2) with the distribution of eq. (4.3.3) so that the same point distribution along the r? coordinates in the initial grid would be retained in the grid obtained using the Adaptive Grid code. Each of the figures used in the examples show only the points in the reduced grid network. 61 50 here because of problems occurring when using the Beam and Warming scheme in conjunction with the adaptive grid generation technique. These problems are discussed more fully in Chapter VII. The first flow problem considered is a the axisymmetric projectile flow problem with an attached sting to eliminate the base flow region. In both the Thin Layer code and the TVD code the boundary conditions are implemented explicitly by updating the flow variables at the boundary after each time step. Referring to Figure 5.2 which shows the basic C- type grid configuration for the axisymmetric projectile flow problem, there are four different sets of conditions applied at the four bounding coordinates. Along the outer boundary, line AB, free stream conditions are assumed. Thus the outer boundary is required to be sufficiently far from the projectile so that this assumption is valid. Along the upstream line of symmetry, line AD, the variables are obtained by requiring symmetry about the x axis. Numerically this is done by setting a second order forward difference approximation to the first derivative along the Â£ coordinate for each variable equal to zero and solving for the values at the points on the axis in terms of known values in the domain. Along the dawmstream boundary, line BC, the variables are obtained by extrapolation along the Â£ coordinate lines. However, if the free stream velocity is subsonic the pressure is fixed. This is done to prevent oscillations in the pressure distribution that would otherwise occur [24]. Along the projectile surface, line CD, the no slip condition holds for viscous flow, thus the velocity components are zero. The density is obtained by zero order extrapolation along the rÂ¡ coordinate lines and the pressure is obtained by satisfying the momentum equation along the surface. 10 It is advantageous to designate the curvilinear coordinates with integer values. Finite difference expressions in the computational plane are simplified since AÂ£ and A>? became unity and the mapping of eqs. (2.1.2) will appear as a two dimensional array with integer subscripts, creating a naturally ordered way to store the grid points. Using the inverse mapping, the location of any coordinate line in the physical domain can be determined simply as the set of points obtained by holding one of the arguments in the mapping constant. The integer indicies i and j are used here to designated each grid point. Thus the 1th point along a Â£ coordinate is and the point along an rÂ¡ coordinate is r? j. For simplicity the notation (xi,j,yi,j) is used to represent the point (x(ii,?j),y(Â£i,ij)) and the total number of points in each direction Â£ and r? are IM and JM, respectively. For the projectile problem considered, the Â£ and rÂ¡ coordinate lines will always run in the streamwise and normal directions as indicated in Figure 2.1. Once a mapping is defined, in a discrete sense, the governing equations can be obtained in the computational plane by transforming them onto the curvilinear coordinate system. Derivatives with respect to the Cartesian coordinates can be expressed using eqs. (2.1.1) and the chain rule as a_ ' ax Â£x Vx ' a ' a? a ay . Â£y r, y - J a . dr> . (2.1.3) where Â£X/ Â£y, r;x, and r?y are the metric coefficients of the transformation. They appear as constants in the transformed equations and depend solely on the mapping of eqs. (2.1.1). Substitution of eqs. 18 'Y7 o 2 2 2 ?xx = (y^x^-2y|yf?x^-y^x)+(y^y^-2y^yr?y^+y^y^) (2.3.3) Similar expressions can be derived for each of the other second derivative terms in eq. (2.2.7), they are listed in Appendix A. All the first derivative terms, of course, can be transformed using eqs. (2.1.4). By substituting these expressions into eqs. (2.2.7) and collecting like terms the equations for adaptive grid generation became, in the computational domain, alxÂ£ Â£ + a2XÂ£7 + a3x^^ + biyÂ£Â£ + b2y^^ + b3y^^ = S1 (2.3.4) clxÂ£ + c2 xÂ£r? + c3Yr)r) + dlY^ + ^Y^r, + *3Yr,r, = S2 in which al = -Yn(a + A'Â£2) 2A 1 y^a/3 a2 = 2yrj (p + A *07) + A'y^ (4/32+J2) a3 = ~YrÂ¡(y + A 172) 2\'y^ap bi = X^ (a + A /32 ) 2A 1 yÂ¡:a/3 b2 = -2x^(0 + A'/37) A 1 x^ (4/32+J2) b3 = xr? 7 + A '72) + 2\'y^p (2.3.5) C1 = y^ (a + A a2 ) + 2 A 1 y^a/3 c2 = ~2y^(/9 + A'a/8) A y,, ( 4/3 2+J2 ) c3 = y? (7 + A 'p2) + 2A 'y^yP dl = xÂ£ (a + A 1 a2) + 2A y^a/3 d2 = 2x^(p2 + \'ap) + A 'y^ (4/32+J2) d3 = -y$(7 + A *^2) + 2\'yr)yp 58 configuration. The difference could also be due to poor resolution in the projectile comer region or an inadequate turbulence model near the flow separation region. 70 Figure 6.3a Grid network upstream of the projectile Figure 6.3b Lack of convergence over secant ogive with Beam and Warming scheme 95 Figure 8.7b Adapted grid network for converged solution for Mach 1.2 Figure 8.8 Setting maximum value for the control function 35 Figure 3.5 Corresponding grid point spacing I Figure 3.6 Initial reduced grid network 8 xml version 1.0 encoding UTF-8 REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd INGEST IEID EJJEMFLZN_42J0HG INGEST_TIME 2014-10-03T21:52:12Z PACKAGE AA00025703_00001 AGREEMENT_INFO ACCOUNT UF PROJECT UFDC FILES 38 The methods used to calculate the control functions from the solution variables is described in Chapter IV. 13 result in a set of elliptic equations which govern the grid. To develop an adaptive grid generation scheme using variational techniques, a measure of each grid characteristic is defined and its integral is taken over the physical domain. If the measures are defined such that they measure a deviation from the desired property then the coordinates Â£ and rÂ¡ that minimize the integral will yield a grid with the desired characteristics. As a minimization problem, the Euler-Lagrange equations can be applied to obtain a set of partial differential equations, the solution of which defines the curvilinear coordinates. The first of three functionals used here is I P vf-.-VÂ£ p dxdy (2.2.1) which measures the adaptation of the grid spacing in the Â£ coordinate direction to the control function P. When P is small, the quantity VÂ£ must also be small in order to minimize the integral; a small value of VÂ£ corresponds to large grid spacing. Consequently, when P is large the grid spacing along the Â£ coordinate will be small. Similarly, the integral I q Vf? Vr> Q dxdy (2.2.2) is a measure of the adaptation of the grid spacing in the rÂ¡ direction to the control function Q. The third functional In = (VÂ£.Vr?)2 dxdy (2.2.3) CHAPTER V FLCW SOLVERS V. 1 The Governing E&uations The governing equations for viscous transonic flow are the Navier- Stokes equations and an equation of state. The Navier Stokes equations, which represent the conservation of mass, momentum and energy form a set of parabolic/hyperbolic partial differential equations when the time dependent terms are retained. To solve these equations numerically on an arbitrary domain the governing equations are scaled and then transformed onto the curvilinear coordinates. For problems in three dimensions, three curvilinear coordinates are required, denoted here as Â£, r) and f. However, for bodies of revolution, such as the axisyrametric projectiles at zero angle of attack considered here, the solution is assumed invariant in the circumferential direction, designated as the f direction, and the equations can be reduced since only two independent spatial variables are required. A further reduction in the equations comes from the thin layer approximation. In the thin layer approximation all the viscous derivatives in the direction along a solid surface (the Â£ direction in the present applications) are dropped from the equations. Due to restrictions on computer storage and CRJ time, it is not usually feasible to resolve the viscous terms in both the directions normal to the surface and along the surface. For high Reynolds number flows, such 44 43 s = s(r?) (4.3.1) After solving for the control function, eq. (2.3.10) becomes 1 (4.3.2) The control function then can be determined from the prescribed distribution s by using central difference approximations to evaluate . As noted previously, the resulting distribution may vary somewhat due to the coupling of the two dimensional adaptive grid equations and effects of orthogonality. This approach, which is used in some of the cases solved here to specify the spacing in the normal coordinate direction, is similar to the technique of Thompson et al. [15] to obtain boundary point distribution in the domain interior. The equation for the grid point spacing used to obtain the arclength distribution is Asj As (1+e)n (4.3.3) where e is determined so that Sp matches the total length of the coordinate line (see reference 16), As* is a specified spacing at the projectile surface, and n is equal to the index j. The actual choices for the functions h and g in eqs. (4.1.4) are discussed in Chapter VIII. 79 variables at the vertices of the triangle as (qiAf + <%2a2 + <33 A3) q = (7.3) (Ai + A2 + A3) where A3, A2 and A3 are the areas of the triangles formed with the grid point in the adapted grid network. The basic algorithm for incorporating the adaptive grid generation technique into both the Thin Layer code and the TVD code consists of the following steps: 1. Input an initial grid network into the flow solver code and begin the time-marching algorithm. 2. Advance the solution a specified number of time steps, and submit the current grid network to the Adaptive Grid code. 3. Obtain an adapted grid network using the Adaptive Grid code. 4. Interpolate the flow variables onto the adapted grid network which new becomes the current grid network. 5. Repeat steps 2 through 5 until the solution converges for steady flew problems or repeat the steps until a specified time is reached for unsteady flow problems. This process which couples the flew solver and the adaptive grid generation technique is referred to here as the self-adaptive computational technique. For the two problems considered here, the axisymmetric transonic flew over a projectile with sting and the transonic projectile base flow problem, only steady flew solutions are sought. Since the adapted grid network depends on the flew variables via the control functions, the convergence of the solution implies also that the grid network ceases to change with time. In all the examples of the following chapters the grid networks of Figures 5.4 and 5.10 are used as A SELF-ADAPTIVE COMPUTATIONAL METHOD FOR TRANSONIC TURBULENT PROJECTILE AERODYNAMICS BY CHRISTOPHER WILLIAM REED A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FIORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FIORIDA 1987 To Shelley and Drew ACKNOWLEDGEMENTS The support of Dr. Chen-Chi Hsu has been most encouraging during the course of this work. His guidance and willingness to provide it are appreciated. In giving me the freedom to experiment, discover, fail and succeed I have learned as much about research as I have about the topic. For this I am grateful. My appreciation is extended as well to the committee members, Dr. Leadon, Dr. Kurzweg, Dr. Mikolaitis and Dr. Wilson. This work was partially supported by the APOSR. Computer time was provided by NASA Ames, NSF and the Pittsburgh Supercomputer Center. iii TABLE OF CONTENTS ACKNOWLEDGEMENTS ABSTRACT vi CHAPTER I INTRODUCTION 1 II THE ADAPTIVE GRID GENERATION EQUATIONS 7 II.1 The Curvilinear Coordinate System 7 11.2 A Variational Approach 12 11.3 Inversion of the Grid Generation Equations 17 III NUMERICAL SOLUTION OF THE ADAPTIVE GRID GENERATION EQUATIONS 22 III.1 The Newton-Raphson Iterative Method 22 111.2 Modified Solution Procedure 26 111.3 The Adaptive Grid Generation Code 37 TV CONTROL FUNCTIONS 39 IV.1 The Basic Form 39 IV. 2 Smoothing 41 TV. 3 Other Control Functions 42 V FIOW SOLVERS 44 V.l The Governing Equations 44 V.2 Solution Algorithms for the Governing Equations ... 46 V.3 Examples on Fixed Grid Networks 51 VI PRELIMINARY TESTS AND MODIFICATIONS 61 VI. 1 The General Procedure 61 VI.2 Orthogonality 62 VI. 3 Curvature 62 VI.4 Internal Grid Boundaries 69 VII THE SELF-ADAPTIVE COMPUTATIONAL METHOD 76 iv VIII RESULTS AND DISCUSSION 81 VIII.1 Choices for the Control Functions 81 VIII.2 Results Using the Self-Adaptive Computational Technique: The Beam and Warming Scheme 84 VIII.3 Results Using the Self-Adaptive Computational Technique: The TVD Scheme 90 VIII.4 Applications to the Projectile Base Flow Problem 97 IX OCNCIUDING REMARKS 102 APPENDIX 104 REFERENCES 105 BIOGRAPHICAL SKETCH 108 V Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy A SELF-ADAPTIVE COMPUTATIONAL METHOD FOR TRANSONIC TURBULENT PROJECTILE AERODYNAMICS By Christopher William Reed August, 1987 Chairman: Chen-Chi Hsu Major Department: Engineering Sciences An adaptive grid generation procedure is developed for viscous flow problems. The equations governing the adaptation are based on a variational statement resulting in a set of elliptic equations in which adaptation can occur independently in each coordinate direction. The method allows for explicit control of adaptation and orthogonality while grid smoothness is implicit in the elliptic equations. The adaptive grid generation equations retain a simple relationship between the control functions and the grid point spacing, are capable of efficiently providing the extremely refined mesh in the boundary layer regions and can maintain a smooth grid network in complex geometries. A self-adaptive computational technique has been developed by coupling the adaptive grid generation procedure with a thin layer Navier-Stokes code and a TVD code. In solving an axisymmetric turbulent vi transonic projectile flow problem, the grid network was adapted to the pressure distribution in the streamwise direction and to the velocity distribution in the direction normal to the projectile surface. Results obtained for Mach 0.91, 0.96 and 1.2 indicate that the procedure can provide good adapted grid networks provided good choices are made for the control functions. vii CHAPTER I INTRODUCTION The study of boundary fitted grid systems has became an important part of computational fluid dynamics and, in fact, their development has been cited as a major pacing item in the progress of computational fluid dynamics [1]. The use of such systems makes the application of finite difference techniques versatile and effective in solving complex fluid dynamic problems in arbitrary domains. Boundary fitted grid generation can be viewed as the development of a general curvilinear coordinate system in which the coordinates coincide with the boundaries on the edges of the physical domain. This feature of the boundary fitted grid systems facilitates the application of finite difference approximations at the boundaries and eliminates the need for interpolation between grid points near the boundary. This is particularly advantageous when the boundaries are highly curved or have slope discontinuities since the use of interpolation near these boundaries may have a significant effect on the solution [2]. In the domain interior, the intersection of each family of coordinate lines denotes each grid point and therefore defines the computational mesh. The distribution of the coordinate lines in the physical domain may 1 2 be concentrated to provide high resolution in certain areas but will always form an equally spaced rectilinear coordinate system in the transformed plane. Thus, once the governing equations are transformed onto the curvilinear coordinates the finite difference algorithms used to solve these equations may be developed on the equally spaced rectilinear coordinate system which is inherent in the transformed plane. Certain characteristics of the general curvilinear coordinate system have been shown to affect the accuracy of the finite difference approximations in the transformed plane [3]. These include orthogonality of intersecting coordinate lines and the smoothness of the grid spacing. It is also known that a refined mesh is needed in regions of large physical gradients to obtain accurate approximations. In fact, poor resolution in these high gradient regions can be detrimental to solution accuracy and to the convergence of finite difference algorithms for fluid flow problems [4]. Thus, the development of a good adapted grid network is important for the successful application of finite difference methods to solve complex fluid dynamic problems. The self-adaptive grid generation technique has became an important area in computational fluid dynamics since it has been shown to provide good grid networks for the complex flow fields occurring in transonic and supersonic flows [5,6,7]. The use of boundary fitted curvilinear coordinate systems with transformed governing equations leads naturally to the concept of solution adaptive grids. As constraints of computer storage and CRJ time prohibit uniform mesh refinement throughout the entire domain, the coordinate spacing in the physical domain is varied to increase resolution in only those regions in which the solution is changing rapidly. For simple flow problems, when the position of 3 important solution gradients is known, good adapted grids can be obtained with conventional techniques. However, in more complex problems the position of these important regions is not known a priori and then the development of a proper adaptive grid network is difficult. Solution adaptive grid generation adresses this problem by continuously updating the grid network as the solution evolves such that the important physical gradients are sufficiently resolved as they develop. It is thus an attempt to efficiently reduce the truncation error by only refining the mesh in those regions where the error may be large. Also, since the grid characteristics of smoothness and orthogonality also affect the truncation error, a good adaptive grid generation technique should include the enhancement of these characteristics as well. The general approach of most adaptive grid generation schemes is based on minimization techniques. A measure of each desired grid characteristic is defined, and the governing equations for the grid are obtained by minimizing the integral of these measures over the domain. In many instances adaptation is required in only one coordinate direction. Ewyer et al. [8], for instance, have developed an equidistribution law to control the grid spacing along one coordinate direction. The spacing was adapted to a flame front in a two-dimensional combustion problem and the improved resolution eliminated spacial oscillations in the flame front. Gnoffo [9] has developed another one dimensional adaptive scheme in which tension springs are assumed to connect the grid points along a coordinate. The spring constants became large when increased resolution was required, and the points clustered in those regions in an attempt to minimize the potential energy of the spring system. The application of this technique to supersonic flow over a Viking Aeroshell successfully 4 resolved a separated shear layer by adapting to the velocity gradient. One dimensional adaptation can be extended to higher dimensions by successive adaptation in each coordinate direction. Nakahashi and Deiwert [7] have used the spring analogy in such an approach and added torsional springs that connect intersecting coordinate lines to control orthogonality. They solved a transonic viscous airfoil problem in which the grid was adapted to the density gradient in the streamwise direction to resolve shocks and adapted to the velocity gradient normal to the airfoil surface to resolve the boundary layer. The one-dimensional approach to multi-dimension adaptation has the advantage of efficiency since the grid can usually be obtained using a marching solution algorithm. Another advantage is the independence of adaptation in each coordinate direction. However, it is difficult to maintain smoothness in such approaches, and highly skewed grids may result [10]. Rai and Anderson [11] have used a different approach for two dimensional adaptation. In their scheme each grid point induces a velocity, either repulsive or attractive, on every other point. By choosing some flow gradient to indicate the need for point clustering, each point with a gradient larger than the average gradient attracts points, and those with values below the average repel points. By summing the induced velocities from each point, the velocity of each point is determined and can be integrated over a time step to determine the new point location. In applying this scheme to the two dimensional linearized viscous Burger's equation, they adapted the grid to the derivative of the dependent variable in each direction and showed that such an approach will increase solution accuracy. 5 Saltzman and Brackbill [5] have developed an adaptive grid generation scheme based on a variational approach. A functional is defined to measure each grid characteristic and the minimization of these functionals results in a set of partial differential equations that govern the grid point spacing. With the proper choice of parameters the equations are elliptic, which guarantees a one-to-one mapping and results in a smooth grid. They solved an inviscid supersonic flow problem by adapting the grid cell size to a function of the pressure gradient to capture shocks. However, it is not possible in this approach to adapt the grid spacing independently in each coordinate direction, a feature that is prohibitive if the spacing in different coordinate directions vary by a large amount. The adaptive grid generation technique presented here is similar in approach to that of Brackbill and Saltzman but is developed for applications to viscous transonic flow problems. The solution of these problems usually contains shock structures aligned with one coordinate direction and boundary shear layers parallel to a streamwise coordinate. Also, the grid point spacing can vary by orders of magnitude along different coordinates and it is thus necessary to adapt the grid independently in each coordinate direction. Therefore, a functional is defined for independent adaptation in each coordinate direction such that the resulting equations are elliptic. The use of elliptic equations in grid generation offer some important advantages that motivate their use. They constitute a boundary value problem and can thus be used in any domain. They guarantee a one to one mapping and also, grid smoothness is inherent in the equations. In addition to the equations governing the grid, there are two other topics that require attention when developing 6 an adaptive grid generation technique. First is the definition of the control functions, which are the part of the grid generation equations that dictate the need for adaptation. They are dependent upon the solution and thus provide the link between the grid generation and the governing equations. Also, a method is needed to account for the effect of the moving grid on the solution. Since the solution is defined at each point, any movement will alter the solution. The proposed adaptive grid generation technique is used to adapt the grid for a transonic axisymmetric projectile flew problem which is solved using both an implicit factorized algorithm and a TVD scheme for the thin layer Navier Stokes equations. The constraint of axisymmetric flow reduces the domain to two independent spatial variables, thus the grid generation equations are developed in two dimensions. The technique, however, can be readily extended to three dimensions. CHAPTER II THE ADAPTIVE GRID GENERATION EQUATIONS II. 1 The Curvilinear Coordinate System The use of boundary fitted curvilinear coordinate systems in solving equations on arbitrarily shaped domains was originally motivated due to problems in applying finite difference approximations at curved boundaries. Because curved boundaries do not in general coincide with the rectilinear Cartesian coordinates, the grid points defined along these coordinates do not usually lie on the bounding curve. It therefore becomes difficult to accurately represent the boundaries [2] and special procedures must be used in those regions [12]. A solution to this problem, which has become an effective and versatile approach, is to map the arbitrarily shaped physical domain (x,y) into a rectangular computational domain (Â£, r?). Such a mapping for a typical projectile flow domain is shown in Figure 2.1. As indicated, the curved Â£ and rÂ¡ coordinate lines coincident with the boundaries in the physical domain map onto the edges of the rectangular computational domain. The boundary grid points, being defined by the coordinate intersections, are coincident with the boundary in the computational plane and thus no special treatment for the boundary conditions is required. Another advantage is that virtually any set of nonuniformally spaced curvilinear coordinates in the physical domain map into an equally spaced regular 7 8 9 mesh in the computational domain. Furthermore, once a successful finite difference algorithm is developed on the computational domain, any physical domain may be mapped into this region and solved with the same numerical algorithm. The task of numerical grid generation is to define the transformation between the Cartesian and computational coordinate systems which can be viewed as the development of a curvilinear coordinate system in the physical domain. The general mapping between the two coordinate systems is = (x,y) (2.1.1) 1 = v (x,y) Owing to the discrete nature of finite difference approximations, it is sufficient to define a finite set of coordinate lines, the intersections of which define the grid points in the physical domain. Thus for each grid point (x,y) in the physical domain, there are two curvilinear coordinates (Â£,??) that intersect at that point and designate its location. It is convenient to view this mapping in the computational domain, where for each pair (|,r?) there corresponds a point (x,y) in the physical domain. This correspondence can be conveniently expressed as the inverse mapping x = x (Â£ rÂ¡) Y = (2.1.2) 10 It is advantageous to designate the curvilinear coordinates with integer values. Finite difference expressions in the computational plane are simplified since AÂ£ and A>? became unity and the mapping of eqs. (2.1.2) will appear as a two dimensional array with integer subscripts, creating a naturally ordered way to store the grid points. Using the inverse mapping, the location of any coordinate line in the physical domain can be determined simply as the set of points obtained by holding one of the arguments in the mapping constant. The integer indicies i and j are used here to designated each grid point. Thus the 1th point along a Â£ coordinate is and the point along an rÂ¡ coordinate is r? j. For simplicity the notation (xi,j,yi,j) is used to represent the point (x(ii,?j),y(Â£i,ij)) and the total number of points in each direction Â£ and r? are IM and JM, respectively. For the projectile problem considered, the Â£ and rÂ¡ coordinate lines will always run in the streamwise and normal directions as indicated in Figure 2.1. Once a mapping is defined, in a discrete sense, the governing equations can be obtained in the computational plane by transforming them onto the curvilinear coordinate system. Derivatives with respect to the Cartesian coordinates can be expressed using eqs. (2.1.1) and the chain rule as a_ ' ax Â£x Vx ' a ' a? a ay . Â£y r, y - J a . dr> . (2.1.3) where Â£X/ Â£y, r;x, and r?y are the metric coefficients of the transformation. They appear as constants in the transformed equations and depend solely on the mapping of eqs. (2.1.1). Substitution of eqs. 11 (2.1.3) into the governing equations makes Â£ and rÂ¡ the independent spatial variables and now a numerical algorithm can be developed using finite difference approximations on the equally spaced computational domain. Since the metric coefficients depend only on the mapping, virtually any domain may be mapped into the computational domain and solved using the numerical algorithm. The metric coefficients can be evaluated for a numerically generated transformation by considering the metrics of the inverse transformation. The relationship between the metric and inverse metrics is, for two dimensions, Â£ x ?y 1 ' Yr, Y ' >1 >1 * J . ~xn . where J is the Jacobian of the transformation, J = x^Yr, ~ xnYt (2.1.5) The derivation of these equations can be found in reference [13]. Since Â£ and rÂ¡ are the independent variables for the inverse metrics, the inverse metrics can be evaluated numerically in the computational domain and subsequently used to determine the metric coefficients. The primary disadvantage of using boundary fitted curvilinear coordinate systems is the appearance of additional truncation error terms in the finite difference approximations on the computational plane. Mastin [3] has shown, for instance, that the truncation error T for a second order finite difference approximation of the first derivative is, 12 in the computational domain T ~ "gX| fxxx2 xÂ£Â£ fxx (2.1.6) The first term is the usual truncation error term containing the grid point spacing x^ and the third derivative of the function. The second term is due to the varied spacing along the Â£ coordinate lines in the physical domain. It is thus possible to increase the truncation error if the coordinates spacing is not smooth. They have also shown that truncation error for a first derivative varies inversely as the sine of the angle between the curvilinear coordinates and therefore a lack of orthogonality can also increase the truncation error. Orthogonality is also important near solid boundaries for turbulent flow problems because turbulence models often are based on information in the normal direction. In the context of a coordinate transformation, the adaptive gridding can be viewed as an attempt to reduce the truncation error by adjusting the transformation metrics. This is accomplished by judiciously defining the curvilinear coordinates such that they vary smoothly, are nearly orthogonal, and concentrate in regions where the solution is changing rapidly. II.2 A Variational Approach The choice of a method to generate the curvilinear coordinates is arbitrary in that it has no dependence on the partial differential equations that are to be solved. Of the many methods available, however, the variational approach is attractive for purposes of adaptive grid generation because it can provide a simple, explicit means of incorporating desired properties into a grid generation scheme and can 13 result in a set of elliptic equations which govern the grid. To develop an adaptive grid generation scheme using variational techniques, a measure of each grid characteristic is defined and its integral is taken over the physical domain. If the measures are defined such that they measure a deviation from the desired property then the coordinates Â£ and rÂ¡ that minimize the integral will yield a grid with the desired characteristics. As a minimization problem, the Euler-Lagrange equations can be applied to obtain a set of partial differential equations, the solution of which defines the curvilinear coordinates. The first of three functionals used here is I P vf-.-VÂ£ p dxdy (2.2.1) which measures the adaptation of the grid spacing in the Â£ coordinate direction to the control function P. When P is small, the quantity VÂ£ must also be small in order to minimize the integral; a small value of VÂ£ corresponds to large grid spacing. Consequently, when P is large the grid spacing along the Â£ coordinate will be small. Similarly, the integral I q Vf? Vr> Q dxdy (2.2.2) is a measure of the adaptation of the grid spacing in the rÂ¡ direction to the control function Q. The third functional In = (VÂ£.Vr?)2 dxdy (2.2.3) 14 measures the orthogonality of the Â£ and rÂ¡ coordinates. As will be shown later, grid smoothness can be controlled through the definitions of the control functions. These three integrals can be combined into one total functional lip + hullL + a (VÂ£ .V7) 2 \ dxdy (2.2.4) p Q J where A is a parameter that weighs the relative importance of adaptation and orthogonality. Before applying the Euler-Lagrange equations for f and rÂ¡, it is beneficial to scale each term in the equations so that they remain of the same order of magnitude. An ordering analysis shows that the first term in eq. (2.2.4) is of order (C2/P), the second is of order (C2/Q) and the orthogonality term is of order (C4/!,2) where C is a computational length scale and L is a physical length scale. It is common practice to use highly stretched grids in transonic flow problems to provide extremely small grid spacing in the viscous sublayer region. The spacing along the normal coordinate can increase by five orders of magnitude, consequently the terms in eq. (2.2.4) may vary in a similar way. It is therefore beneficial to scale the terms locally rather than globally since a global scale cannot reflect the size of the term throughout the domain. The local scales Ap, Aq, and A0 used here are Ap=pL Aq=QL ao=jL (2.2.5) where Pp and Qp are the local values of the control functions and Jl is the local Jacobian and is of order (L2/C2). Although they vary throughout the domain they are considered constants during the application of the 15 Euler-Lagrange equations. This violates the variational principal, but it has been shown [13] that the use of local scaling will produce better grids when there are large changes in the grid spacing. After substituting in the scales, the total functional Ip becomes It= [ { Pl V^:VÂ£ + QlV>? Vr? + ajL (v^.Vr?)2 }dxdy (2.2.6) After applying the Euler-Lagrange equations for Â£ and rÂ¡ to the total functional I c + c + yp-.ve xx + + ^ (*?xÂ£xx+2r7x7yÂ£xy+riyyy + + (v xÂ£ y+r) yÂ£ x) 7 xy+ (Â£ x7 x+2Â£y7y) 7yy) (2.2.7) 7XX + 7yy + + (?x7x+2Â£y7y) Â£xy + * ( (2Â£x7x+Â£y7y) ?xx+ (rlx^y+riy^x) Â£xy + 2 p ^x7xx+2^x^y7xy + Cy7yy)=0 The subscript L has been dropped from the local scales in the differential equations since they are no longer needed to indicate their assumed invariance. The Â£ and rÂ¡ coordinates that satisfy these equations will yield a grid with the desired properties defined in the functionals of eqs. (2.2.1), (2.2.2) and (2.2.3). It should be noted that there was no explicit functional to maintain grid smoothness. There are two types of smoothness that should be differentiated. One type, which can be controlled by the control functions, is smoothness of the grid point 16 spacing along a coordinate line. If the control functions change rapidly along a coordinate, the grid spacing will consequently change rapidly. This smoothing thus can be governed through the definition of the control functions. Another type of smoothness, which is inherent in the elliptic equations, is that between the grid point spacing of adjacent coordinates. If the grid point distribution along two adjacent coordinates varies, the elliptic equations will tend to dampen the variation. Equations (2.2.7) govern the grid coordinates Â£ and rÂ¡ and will be elliptic as long as the parameter A is not too large. They constitute a boundary value problem and it is therefore necessary to define the coordinates along the boundaries. It is also necessary to adapt the boundary points in a consistent manner with the interior points so that the boundary points remain aligned with the interior points. For this purpose, a one dimensional functional analogous to equation (2.2.1) or (2.2.2) is defined to measure adaptation along a boundary coordinate, (2.2.8) where s is the arclength along the boundary coordinate. There is no consideration for orthogonality and scaling is unnecessary. Applying the one dimensional Euler-Lagrange equation yields a second order differential equation, Â£ss Â£ s Ps/P0 (2.2.9) Note that a similar equation can be obtained for boundary adaptation along a rÂ¡ coordinate by replacing f in eq. (2.2.9) with r? and P with Q. This one dimensional equation, together with eqs. (2.2.7), are the 17 equations that govern the grid coordinates in the physical domain. II.3 Inversion of the Grid Generation Equations In order to solve the grid generation equations, they are transformed onto the curvilinear coordinate system. This is done by inverting eqs. (2.2.7) and (2.2.9) to make x, y, and s the dependent variables. The curvilinear coordinates then become the independent variables and the grid generation equations are solved on the computational domain which is equivalent to defining the inverse mapping of eqs. (2.1.2). All the advantages of using the computational domain to solve the governing equations discussed previously will be realized here as well. There are two approaches to inverting the equations. One is to transform the functional Ip onto the computational domain and then apply the Euler-Lagrange equations for x and y. Another approach, used here is to derive the corresponding expressions on the computational domain for each term in eqs. (2.2.7) and (2.2.9) and then to substitute them into the equations. The first term in eq. (2.2.7) can be written using eq. (2.1.3) as Â£xx = (Â£x + 7X ~ ) Â£x (2.3.1) dr) Using eq. (2.1.4) to transform the metrics this expression becomes Yr, d Yz d Yv Â£ xx -( ) (2.3.2) J J dr) J and after differentiating, the relationship is found to be 18 'Y7 o 2 2 2 ?xx = (y^x^-2y|yf?x^-y^x)+(y^y^-2y^yr?y^+y^y^) (2.3.3) Similar expressions can be derived for each of the other second derivative terms in eq. (2.2.7), they are listed in Appendix A. All the first derivative terms, of course, can be transformed using eqs. (2.1.4). By substituting these expressions into eqs. (2.2.7) and collecting like terms the equations for adaptive grid generation became, in the computational domain, alxÂ£ Â£ + a2XÂ£7 + a3x^^ + biyÂ£Â£ + b2y^^ + b3y^^ = S1 (2.3.4) clxÂ£ + c2 xÂ£r? + c3Yr)r) + dlY^ + ^Y^r, + *3Yr,r, = S2 in which al = -Yn(a + A'Â£2) 2A 1 y^a/3 a2 = 2yrj (p + A *07) + A'y^ (4/32+J2) a3 = ~YrÂ¡(y + A 172) 2\'y^ap bi = X^ (a + A /32 ) 2A 1 yÂ¡:a/3 b2 = -2x^(0 + A'/37) A 1 x^ (4/32+J2) b3 = xr? 7 + A '72) + 2\'y^p (2.3.5) C1 = y^ (a + A a2 ) + 2 A 1 y^a/3 c2 = ~2y^(/9 + A'a/8) A y,, ( 4/3 2+J2 ) c3 = y? (7 + A 'p2) + 2A 'y^yP dl = xÂ£ (a + A 1 a2) + 2A y^a/3 d2 = 2x^(p2 + \'ap) + A 'y^ (4/32+J2) d3 = -y$(7 + A *^2) + 2\'yr)yp 19 = *r,2 + Yb2 7 = XÂ£2 + y^2 P = + xrjYr, A' = A/J and Si = P^Ja P PnJ/9 s2 = Q Qt7J7 + Q (2.3.6) It is important to note that in inverting the equations the derivative of P in the r? direction appears even though P was intended to control the spacing in the Â£ direction. This occurs because the contravariant base vector VÂ£ actually indicates the spacing in the direction normal to lines of constant Â£ rather than in the Â£ direction. If the grid is orthogonal, the two directions, VrÂ¡ and VÂ£, are equivalent and the r? derivative will not have any influence since p will vanish. It has been found in numerical experiments that eliminating this term from the equations has a negligible effect and consequently it is dropped from the equations. The source terms S1 and S2 in eqs. (2.3.5) then become P^Ja Qt7J7 (2.3.7) S2 = Q 20 The one dimensional equation for the boundary adaptation, eq. (2.2.9), is also transformed in a similar manner. The second derivative term can be shown, using the chain rule, to transform as sÂ£ Â£ ss = T (2.3.8) sr Substituting eq. (2.3.7) into eq. (2.2.9) yields the transformed one dimensional equation for adaptation along Â£ boundaries + sÂ£ PÂ£/P =0 (2.3.9) and for adaptation along boundaries coincident with rÂ¡ coordinates s^ + s^Q^/Q = 0 (2.3.10) These equations can be solved directly for the grid point spacing along a coordinate. For instance, integrating eq (2.3.10) once yields the grid point spacing s r, = c/q (2.3.11) where c is a constant of integration and can be shown to be C = (2.3.12) f da. J Q where Sp> is the total arclength along the coordinate in the physical domain. However, since the solution in the interior will require an iterative algorithm, the differential form, eqs. (2.3.9) and (2.3.10) are 21 retained and solved iteratively so that the boundary points remain consistent with the interior points as the grid evolves. CHAPTER III NUMERICAL SOLUTION OF THE ADAPTIVE GRID GENERATION EQUATIONS III.l The Newton-Raphson Iterative Method Equations (2.3.3), (2.3.9) and (2.3.10) constitute a set of nonlinear elliptic partial differential equations. In order to solve these equations numerically they are approximated with second order central difference expressions which results in a set of coupled algebraic equations for each grid point. The finite difference approximation for the equations at grid point, (Xj^j, Yifj) become al(xi+l,j 2xi,j + xi-l,j) + a3(xi,j+1 "2xi,j + xi,j-l) + bl(Yi+l,j -2yij +Yi-l,j) + b3(Yi,j+1 -2Yi,j + 2Yi,j-l)-Fi=0 (3.1.1a) cl(xi+l,j "2xi,j + xi-l,j) + c3(xi#j+1 -2xi/j +Xifj_i) + dl(Yi+i,j -2Yi,j +Yi-l,j) + d3 (Yi, j+1 -2Yi,j +Yi, j-l) -F2=:0 (3.1.1b) P|Ja Fi = a2Ti b2T2 (3.1.1c) Qr;J7 F2 = a2Tx d2T2 22 23 T1 = (xi+l,j+l + xi-l,j-l xi+l,j-l ~ xi-l,j+l)/4 (3.1.Id) T2 = (Yi+1,j+1 + Yi-1,j-1 Yi+1, j-1 Yi-1,j+l)/4 where the source terms and the coefficients all contain first derivatives which vhen approximated with central differences do not contain the point xi,j' Yi,j explicitly. This set of equations is solved using a Newton- Raphson iterative scheme. The initial guess for the grid point locations will not satisfy the finite difference representation of the equations and therefore a non-zero residue will exist (Rl)l,j = eq (3.1.1a) (3.1.2) (R2)i,j = eq (3.1.1b) where i. indicates the iteration number. An estimate of the residue at the next iteration can be obtained by expanding the residue in a Taylor series about the current residue and retaining only the linear terms, 3(Rl)i,j a^i,j l a 0 4-1 o a cd- \ . (3.1.3) Â¡L d (r2)i,j aYi,j where ax and Ay are the changes in the grid point location 24 + 1 Z Ax/j xi,rxi,j (3.1.4) 2 + 1 2 AYi/j Yi, jYi,j The derivatives of the residues with respect to the variables and Yi j can be obtained by differentiating eqs. (3.1.1) and are 3(Rl)i,j = -2(ai+a3) , 3xi, j 3(Rj)i,j an,j = -2(b^+b3) (3.1.5) 3(R2),j = -2(C1+C3), 3xi, j 3(R2>,j 3Yi,j -2(d^+d3) An approximate value of AXj^j and AYi,j can be obtained if we require the residue at the next iteration to vanish 2 1 r 3(ri),j 3(Ri),j i/D 3xi, j aYi,j Axi, j 2 3(R2) i, j 3(R2)1,j . Ayi,j . i/j J 3xi/ j aVi, j (3.1.6) Substituting eqs. (3.1.5) into eqs. (3.1.6) and solving for Axj^j and Ay-L, j yields where r Axi,j 1 l *Yi,j J 2 ' (di+d3) -(bi+b3) D . -(C1+C3) (ai+a3) , . (R2)i/j . (3.1.7) D 4((aj+a3)(d^+d3)-(c^+c3)(d^+d3)) (3.1.8) 25 Knowing Axj_j and Ay^j, the new point location can be calculated by . Â£+1 SL+1 solving eqs. (3.1.4) for xj^j and yj^j. This same procedure can be applied to the one dimensional equations for adaptation, eqs. (2.3.9) and (2.3.10). The residue for eq. (2.3.9) is (Rs)i = si+l 2si + si-l + sÂ£pÂ£/p (3.1.9) and the equation for the new point location becomes s{+1=sf+As{ (3.1.10) As{=(Rs)i/2 (3.1.11) The numerical algorithm consists of updating the value of each point on the boundary with eq. (3.1.11) and of sweeping through the domain and updating each interior point using eq. (3.1.7). The criterion used here to indicate convergence is Rmax ^ s (3.1.12) where F^ax is the maximum residue in the interior weighted by the Jacobian Axl,j + Ayf j Rmax= max ( )1/2 (3.1.13) Currently, the criterion 6=0.005 is used. 26 In the original iterative procedure, the most recently updated values for each grid point location were used in updating the remaining points and thus the procedure resembles a Gauss-Seidel iteration. This approach, however, is not suitable for vector processing and thus cannot benefit from the increased processing speed. Therefore, another approach is used in which all the grid points along a coordinate line are updated in a Jacobi iterative fashion. The new grid point values along that coordinate are then used in updating the points along the next coordinate line. Thus the technique can be viewed as a half Jacobi, half Gauss-Seidel iteration. This technique, which is suitable for vector processing, reduced the computer processing time for each iteration by 40 percent. Furthermore, only a few more iterations were required for the same convergence criterion in comparison to the Gauss-Seidel iteration, thus the effective speed up in processing time is approximately 40 percent. Ill.2 Modified Solution Procedure As pointed out earlier, the grids used for viscous transonic flow problems contain highly refined grid spacing in the direction normal to the surface (rÂ¡ direction) in order to resolve the viscous sublayer. This spacing results in grid cells with aspect ratios, up to 105 in the present applications, near solid boundaries. Iterative solution procedures for elliptic equations become inefficient for such large aspect ratios since the motion in one direction will be limited by the small spacing in the other direction. The effect of the aspect ratio on an iterative process can be ascertained by considering the simple grid cell shown in Figure 3.1 which for convenience is rectangular and aligned with the Cartesian Figure 3.1 High aspect ratio cell Figure 3.2 Initial grid network 28 coordinates. For this example the aspect ratio AR is equal to (1/b). Assume for purposes of demonstration that the two terms (Qv/Q) and (P^/P) are non-zero at the center point and will therefore cause the point to move. We can obtain the changes in the point location by evaluating eqs. (3.1.7) for this simple case. Many terms for this example are zero xr?=o yz=0 y?7=0 y^=0 x*e= y^=0 xr,r,=0 Yr,r,=0 The coefficients became a3=b3 b]_=0 c3=0 d3=b2 a2=0 b2=0 c2=0 d2=0 a3=b b3=0 c3=0 d3=l (3.2.1) (3.2.2) and upon substitution into eqs. (3.1.7) yield Axi, j ( JL 2P l+b AYi,j = 2Q l+b (3.2.3) Writing these in terms of the aspect ratio AR, its effect on the grid point motion becomes clear. If we compare the solution for a square cell , i.e. AR=1, 29 Axi,j = 2P 2 ^r? 1 Ayi^ = V 2Q with that for a cell of large aspect ratio AR 1, mj.2 (-L-, L'-1 2P 1+AR (3.2.4) _1 AR . Ayi -i ( ) -L,J 2Q 1+AR (3.2.5) it becomes evident that the reduction in the grid point motion in the x direction for large AR is proportional to (1/AR). For aspect ratios of the order (105) this result becomes prohibitive, rendering the iterative procedure extremely inefficient. An example has been formulated to demonstrate this deficiency. In an initial grid network, as shewn in Figure 3.2, the points along the Â£ coordinates are clustered in one region and the points along each rÂ¡ coordinate are clustered near the bottom surface to form the highly stretched spacing typical of grid networks used in transonic flow calculations. The spacing at the surface is approximately 105 times that in the upper region of the grid. In the adaptive grid generation equations, the control function P is set equal to a constant so that the points will redistribute themselves with equal spacing along the Â£ coordinates. Using a technique described in Chapter IV, the control function Q is defined so that the highly stretched spacing along each rÂ¡ 30 coordinate will be retained in the solution to the adaptive grid generation equations. Using the Newton-Raphson iterative solution procedure, the solution is advanced 20 iterations which results in the grid shown in Figure 3.3. As can be seen, the point spacing along the Â£ coordinates is approximately equal for the coordinates sufficiently above the surface. However, in the region of the high aspect ratio cells, the movement of the points in the Â£ direction is retarded by the extremely small spacing in the x] direction. A modification has been made to eq. (3.1.11) which updates the points along the bounding Â£ coordinate in formulating this example. The movement of points lasing the one dimensional equation is independent of the interior points and thus is not restricted by the high aspect ratio cells. This creates an inconsistency at the boundary since the points adjacent to the boundary are constrained by the high aspect ratio cells. Equation (3.1.11) therefore has been modified as i a i ASi = (Rsi)i ( ) (3.2.6) 1+AR. l where AR* is the aspect ratio of the cell adjacent to the 1th point along the boundary Â£ coordinate. This change only effects the rate of convergence of the iterative solution procedure for the boundary points so that their movement will remain consistent with that of the interior points adjacent to the boundary during the convergence process. In order to alleviate the deficiency of the above solution procedure, a technique has been developed which is based on solving for a reduced grid in which many of the points in the r? coordinate direction have been 31 Figure 3.3 Grid network after 20 iterations f Figure 3.4 Reduced aspect ratio cell 32 removed to decrease the maximum cell aspect ratio. Figure 3.4 shows a diagram of a typical grid for viscous flow problems containing large aspect ratio cells near a solid boundary. A reduced grid can be formed by removing all the Â£ coordinate lines denoted with the dashed lines, resulting in a grid with fewer points along the r? coordinate lines which are designated by the intersection of the solid lines. Currently, enough points are removed along the rÂ¡ coordinate lines such that the minimum spacing between Â£ coordinate lines is greater than 0.04. The typical spacing along the boundary is of order (0.1) which results in cell aspect ratios of order (1). Upon convergence of the solution for the reduced grid, all the points removed are reinserted along coordinate lines using the same one dimensional equation that controls the spacing along the boundaries. It is necessary in this process to adjust the control function Q at the remaining points to account for the spacing requirements of the points removed. A simple procedure for this purpose can be derived by using the equation for one dimensional adaptation, eq. (2.3.10). The procedure will also be adequate for the two dimensional adaptive grid generation equations. If we write explicitly the finite difference approximation for eq (2.3.10), sj+l-sj-l Qj+lQj-l sj+l-2sj + sj_! + ( )( )/Qj=0 (3.2.7) we can regroup the terms in the form Qj+lQj-l Qj + ( ) (sj+l-sj)=(Sj-Sj_1) ( ) (3.2.8) Qj+l"Qj-l Qj ( ) 33 which, after using forward difference expressions to indicate the spacing between two adjacent points can be written ASj = AsjCj (3.2.9) C j + Qj+iQj-i Qj+i~Qj-i 2 (3.2.10) Now, as shown in Figure 3.5, if the ith point is removed it will be necessary to adjust the control function such that the new relationship is satisfied Asj + = ASj_! (3.2.11) where Cj_x is the modified control function and ASj_x is the increased grid point spacing ASj A Sj 4- ASj (3.2.12) The correct form of Cj_x can be obtained in the following derivation "k k Asj + 1 = ASj_x Cj _x ASj Cj + x (ASj_x+Asj) Cj_i Asj 1 Cj Cj+i = Asj _x (1+Cj) Cj_x C1C1 + 1 * = C* 1+C-i J 1 (3.2.13) 34 Once Cj_i is calculated losing eq. (3.2.13), the modified value of Qj-i * can be obtained by solving eq. (3.2.10) for Qj-i 1+C1-1 Q*-l = ( ) (Qj + i-Qj-2) (3.2.14) 1-Cj_i This procedure is applied for each point along an coordinate that is removed. It should be noted that this procedure does not affect the final grid, the points remaining in the reduced grid will obtain the same location as if all the points were used in the iterative solution procedure. The removed points are, once the solution for the reduced grid is obtained, inserted back along the r? coordinates lines. This is done by approximating the r? coordinate with straight line segments between the points of the reduced grid and then reinserting the removed points along each segment with spacing dictated by eqs. (2.3.10) using the original control function Q. If necessary, the use of a higher order approximation to the r? coordinate lines, such as a cubic spline, can be used to obtain a better representation of the rÂ¡ coordinate lines. In using this procedure with the two dimensional grid generation equations, the solution may vary somewhat since the equations are coupled. Also, the effects of orthogonality may alter the results in some areas. In order to show the advantage of this procedure, the initial grid network of the previous example is advanced 20 iterations using the reduced grid technique. Figure 3.6 shows the initial reduced grid in which the points responsible for the high aspect ratio cells cure removed. After modifying the control function Q, the solution to the adaptive grid generation equations is advanced 20 iterations resulting in the grid shown in Figure 3.7. As can be seen, the points along all Â£ coordinates 35 Figure 3.5 Corresponding grid point spacing I Figure 3.6 Initial reduced grid network 36 Figure 3.7 Reduced grid network after 20 iterations Figure 3.8 Final grid network after reinserting grid points 37 are approximately equally spaced and the spacing along the r? coordinate lines in the initial reduced grid is retained which indicates that the modified control function is correct. Upon inserting the removed points along the rÂ¡ coordinates using eq. (2.3.12), the final grid network, Figure 3.8, is obtained, which is the same grid that would have been obtained in Figure 3.3 after a large number of iterations. Ill.3 The Adaptive Grid Generation Code A computer code base on the iterative solution procedure described in the preceding section has been developed and is referred to here as the Adaptive Grid code. As the solution process is iterative, an initial guess for the grid point locations is required and the flow variables for which the grid network is to be adapted cure needed for calculation of the control functions. These are considered as input to the Adaptive Grid code. The initial guess for the grid network can be a grid network obtained by any conventional grid generation technique or a previously adapted grid network. Values for the flow variables at the grid points are obtained using the flow solvers which are described in Chapter V. The sequence of calculations in the Adaptive Grid code proceeds as follows: 1. Input an initial grid network and the flow variables at each grid point. 2. Calculate the control functions P and Q. 3. Remove points to form the reduced grid network and calculate the modified control function Q*. 4. Solve the grid generation equations for the reduced grid network using P and Q*. 5. Reinsert the removed points along the coordinate lines using the control function Q which then results in the adapted grid network. 38 The methods used to calculate the control functions from the solution variables is described in Chapter IV. CHAPTER IV CONTROL FUNCTIONS IV. 1 The Basic Form It is the role of the control functions in the equations governing grid adaptation to dictate the need for increased resolution and, as indicated in section II. 1, they should be judicially chosen so that the resulting grid resolution reduces what would otherwise be large truncation error. For instance, Pierson and Kutler [14] chose for the control function, when solving a one dimensional inviscid Burger's equation, the third derivative of the dependent variable which is also the leading truncation error term in their finite difference approximation to Burger's equation. Their one dimensional adaptive grid scheme was thus an attempt to minimize a measure of the truncation error over the domain. In more complex problems, however, with more dependent variables and in higher dimensions, the truncation error expressions are complex and cumbersome to calculate. Therefore, in practice, an understanding of the physical problem and experience guide the choices for the control functions. The general control function considered here is P = 1 + 7pf Q = 1 + 7qg (4.1.1) 39 40 where 7p and 7q are parameters and f and g are some derivatives of the dependent variables scaled to range between zero and one. In the following developments, the terms Wf7 and h will be used to indicate either P,7p and f or Q,7q and g. By evaluating numerically the solution to the one dimensional equation for adaptation, the following expression equating the control function to the grid point spacing between points (i) and (i+1) along a coordinate line is obtained, ASi Sijr 1 l+7hj S j 1 l+7hj (4.1.2) After evaluating this expression for the maximum and minimum spacing corresponding to h = 1 and h = 0 respectively, and forming their ratio (Asi) max = 1 + 7 (4.1.3) (Asi)mim it can be seen that the specification of 7 dictates the ratio of maximum and minimum spacing along a coordinate line. The parameter 7 also influences the smoothness of the resulting grid point distribution. If it is chosen as zero, equal spacing as indicated by eq. (4.1.2) would result; as 7 is increased, the spacing must change more rapidly to obtain the varied spacing. Following the work of Nakahashi and Deiwert [7], another parameter a has been introduced to obtain more control over the spacing. The general control function now becomes W 1 + 7ha (4.1.4) 41 By prescribing the minimum spacing and 7, eq (4.1.2) can be written as (Asi)min ? (i+7h? ^ 1+7 = 0 (4.1.5) in which a is the unknown quantity and can be determined by a root finding technique. Currently, the Newton-Raphson finding method is employed. Thus, both the minimum and maximum grid point spacing along a coordinate can be specified by specifying the two quantities, (As^J^^ and 7. It should be noted that the prescribed values may not be realized in practice due to the coupling effects of the two dimensional adaptive grid generation equations. Also, the constraint of orthogonality may alter the results in some areas. IV.2 Smoothing The choices for the functions f and g will of course depend on the problem being solved. However, the functions will in general, contain derivatives of the solution which are approximated numerically and it is therefore useful to smooth the control functions to eliminate any spurious data. Also, since smoothing is most effective when the function changes most rapidly, it will help reduce any rapid variations in the control function and subsequently smooth the grid point spacing. The method of smoothing used here, which is similar to that used in reference [5] is 42 2 i=2,IM-1 j=l,JM (4.2.1) 2 i) i=l,IM j=2,JM-1 where the value of w used here is 0.4. The control function at a point is replaced by a weighted average of itself and the two adjacent points along a coordinate line. Currently, two consecutive sweeps of this algorithm are made over the domain for each of the control functions, P and Q. In this algorithm there is no accounting for the varying spacing between points and thus this smoothing is in the computational plane. IV.3 Other Control Functions In many instances a certain point distribution in one coordinate direction is known to yield good results and it would be advantageous to maintain this prescribed distribution in that direction while allowing adaptation to the solution in the other directions. For instance, in some of the viscous transonic flow problems solved here, it is known that a distribution in the direction normal to the surface based on an exponential function will yield good results provided enough points are used. A procedure to incorporate specified point distributions into the adaptive grid scheme is easily developed by using eq. (2.3.10). Consider some given grid point distribution s in the rÂ¡ coordinate direction which is to be reproduced by the adaptive grid generation equations 43 s = s(r?) (4.3.1) After solving for the control function, eq. (2.3.10) becomes 1 (4.3.2) The control function then can be determined from the prescribed distribution s by using central difference approximations to evaluate . As noted previously, the resulting distribution may vary somewhat due to the coupling of the two dimensional adaptive grid equations and effects of orthogonality. This approach, which is used in some of the cases solved here to specify the spacing in the normal coordinate direction, is similar to the technique of Thompson et al. [15] to obtain boundary point distribution in the domain interior. The equation for the grid point spacing used to obtain the arclength distribution is Asj As (1+e)n (4.3.3) where e is determined so that Sp matches the total length of the coordinate line (see reference 16), As* is a specified spacing at the projectile surface, and n is equal to the index j. The actual choices for the functions h and g in eqs. (4.1.4) are discussed in Chapter VIII. CHAPTER V FLCW SOLVERS V. 1 The Governing E&uations The governing equations for viscous transonic flow are the Navier- Stokes equations and an equation of state. The Navier Stokes equations, which represent the conservation of mass, momentum and energy form a set of parabolic/hyperbolic partial differential equations when the time dependent terms are retained. To solve these equations numerically on an arbitrary domain the governing equations are scaled and then transformed onto the curvilinear coordinates. For problems in three dimensions, three curvilinear coordinates are required, denoted here as Â£, r) and f. However, for bodies of revolution, such as the axisyrametric projectiles at zero angle of attack considered here, the solution is assumed invariant in the circumferential direction, designated as the f direction, and the equations can be reduced since only two independent spatial variables are required. A further reduction in the equations comes from the thin layer approximation. In the thin layer approximation all the viscous derivatives in the direction along a solid surface (the Â£ direction in the present applications) are dropped from the equations. Due to restrictions on computer storage and CRJ time, it is not usually feasible to resolve the viscous terms in both the directions normal to the surface and along the surface. For high Reynolds number flows, such 44 45 as the transonic flow considered here, large velocity gradients are expected to occur normal to the surface and most of the grid points are used to resolve these gradients. As larger spacing is then used along the surface and the viscous derivatives in that direction are expected to be small in comparison, they are dropped from the equations. The conservative form of the transformed thin layer governing equations for axisymmetric flow, written in the compact vector form are atq + dtk + dÂ£ + b = (s.i.i) where p P U A : pv , E=J* p w . C , P U pV puU+?xp A X pUV+r?xp pvU+^yp , F=J* pW+flyP pwU+ezp pWV+r7Zp (e+p)U j . (e+p)V 6=J* 0 0 -pw2y p/y pW(y^U+y^V) 0 (5.1.2) A s = J* 0 uRu^ + (u/3)x uRv,, + (u/3)T*y }, uRw^ + (u/3)Tr)z A(%m (u2+v2+w2 ) ^uPf1 (7-I) _1 (a2) ^ + (u/3) (r?xU+r?2yv+r?2xw)T and R = ^x2 + Vy2 + 1z2 T + nyVr, + *?xw>7 46 Here, p is the density, p is the pressure and E is the total specific energy. The velocities u,v and w correspond to the x, y and z directions, respectively, as indicated in Figure 5.1 and U, V, and W are the contravariant velocity components corresponding to the Â£, rÂ¡ and f directions respectively. The coefficient of viscosity is p., Pr is the Prandtl number, Re is the Reynolds number, 7 is the ratio of specific heats and a is the local speed of sound. The details of the transformation and application of the thin layer approximation to obtain eqs. (5.1.1) are available in reference [13]. The assumption of a calorically and thermally perfect gas was made in closing the system of equations. The coefficient of viscosity pi is considered variable and it is also necessary to account for the effects of turbulence for the present applications. Under the assumption of f invariance, only a two dimensional grid network containing Â£ and r? coordinates in a plane through the axis of symmetry will be required. It should be noted, however, that the metric coefficients represent the three dimensional coordinate transformation and are therefore different from the expressions of eqs. (2.1.4) and the Jacobian J* now represents the grid cell volume. V.2 Solution Algorithms for the Governing Equations The implicit spatially factored numerical scheme of Beam and Warming [17] is used to solve the transformed governing equations. The scheme is implemented for three dimensional problems in a code due to Steger and Pulliam [18]; the code used here is a modified version which was developed by Nietubicz et al. [19] for the efficient calculation of axisymmetric projectile problems and is referred to here as the Thin 47 Figure 5.1 Coordinate systems for projectile configuration B Figure 5.2 Boundary configuration for projectile with sting 48 Layer code. The code employs Sutherland's law [20] to relate the viscosity and heat conductivity to the temperature and the effects of turbulence are modeled using the eddy viscosity model of Baldwin and Lomax [21]. The numerical algorithm is a time-marching scheme in which steady state solutions are obtained as the time asymptotic solution of the unsteady equations. Spatial derivatives are approximated with central differences and a forward difference approximation is used for the unsteady terms. The time marching scheme of the finite difference approximation to the governing equations are written in delta form as (I + h S^k ciDi) (I + hfi^ hR^-jfoj1 jD2 ) Aq = h (d^fe + 6n$ ~ Ri1 S fj 6) (5.2.1) e Ed3 where A and B are the Jacobian matrices a = a aÂ£ B = (5.2.2) aq and M is obtained as the time linearization of S. Here, 6 is the central difference operator. Details of this scheme can be found in Warming and Beam's [17] work, as well as in that of others. The terms D3, and D2 are implicit dissipation operators and D3 is an explicit dissipation term Di = J-1 VÂ£ J D2 = J_1 A J D3 = J"1 ((A^ Ve)2 + Ay V,)2) J Â§ (5.2.3) 49 where D and V are forward and backward difference operators, respectively. Since the scheme is not naturally dissipative these terms are added to the finite difference scheme to prevent nonlinear instabilities and high frequency spatial oscillations. The constant coefficients ej and eg used here are 4At and 2At respectively which are consistent with values used in other applications of this code [19]. In general, these coefficients should be chosen large enough to dampen the instabilities but kept small enough to prevent any degradation of solution accuracy. In Chapter VII, seme examples using the Thin Layer code with the Adaptive Grid code shew that there are problems in obtaining a converged solution for the thin layer equations when the grid is adapted. Therefore another numerical scheme for solving the thin layer equations is also used. This scheme is the TVD scheme developed by Yee [22] for the multidimensional Navier-Stokes equations. The theory and development of TVD schemes is beyond the scope of current discussions. However, for purposes here it is sufficient to view the TVD scheme as an improvement in the artificial damping terms in the Beam and Warming scheme. In fact, the TVD scheme developed by Yee for multidimensions can be implemented into the existing Thin Layer code based on the Beam and Wanning scheme by modifying only the added dissipation terms. The original Thin-Layer code has been modified to provide an option for the TVD scheme [23] and is referred to here as the TVD code. These 'sophisticated' dissipation terms switch the scheme from second order to first order accuracy at points of extrema and result in a more dissipative scheme in those regions. It should be pointed out that the primary purpose of the current research is to develop an adaptive grid generation technique. The TVD scheme is used 50 here because of problems occurring when using the Beam and Warming scheme in conjunction with the adaptive grid generation technique. These problems are discussed more fully in Chapter VII. The first flow problem considered is a the axisymmetric projectile flow problem with an attached sting to eliminate the base flow region. In both the Thin Layer code and the TVD code the boundary conditions are implemented explicitly by updating the flow variables at the boundary after each time step. Referring to Figure 5.2 which shows the basic C- type grid configuration for the axisymmetric projectile flow problem, there are four different sets of conditions applied at the four bounding coordinates. Along the outer boundary, line AB, free stream conditions are assumed. Thus the outer boundary is required to be sufficiently far from the projectile so that this assumption is valid. Along the upstream line of symmetry, line AD, the variables are obtained by requiring symmetry about the x axis. Numerically this is done by setting a second order forward difference approximation to the first derivative along the Â£ coordinate for each variable equal to zero and solving for the values at the points on the axis in terms of known values in the domain. Along the dawmstream boundary, line BC, the variables are obtained by extrapolation along the Â£ coordinate lines. However, if the free stream velocity is subsonic the pressure is fixed. This is done to prevent oscillations in the pressure distribution that would otherwise occur [24]. Along the projectile surface, line CD, the no slip condition holds for viscous flow, thus the velocity components are zero. The density is obtained by zero order extrapolation along the rÂ¡ coordinate lines and the pressure is obtained by satisfying the momentum equation along the surface. 51 Another problem considered is the projectile base flew problem. In Figure 5.3, which shows the basic O-type grid configuration, the Â£ coordinate lines bend around the sharp comer at the base and end on a downstream axis of symmetry which is now an rÂ¡ coordinate. The existing codes are used for this problem by simply modifying same of the boundary conditions. The boundary conditions along the upstream axis of symmetry, line AD, and the projectile surface, line CD, are the same as those used in the projectile problem with sting. For the downstream axis of symmetry, line BC, the same conditions along the upstream axis are used; however, a backward finite difference formula is required. The freestream conditions are applied along the outer boundary along the line AA1 and the downstream boundary conditions for the projectile with sting problem are applied along the portion of the outer boundary between A'B, except that the extrapolation is done along the rÂ¡ coordinate lines. V.3 Examples on Fixed Grid Networks In order to show the general characteristics of the Beam and Warming and TVD schemes, and for purposes of comparison with the solutions obtained with the adaptive grid generation technique, the two schemes were used to solve the projectile flow problem with sting on a fixed grid network. The fixed grid network was obtained using a method developed by Steger et al. [16] which is an effective technique for external flow problems. Two equations, one enforcing orthogonality and one specifying the grid cell volume, form a set of hyperbolic partial differential equations which are solved numerically to generate the grid network. In the marching solution algorithm, the initial grid points are specified along the projectile surface. Ihe solution is marched outward in the v 52 Figure 5.3 Boundary configuration for projectile base flow problem SOCBT Figure 5.4 SOCBT projectile configuration 53 coordinate direction, defining each successive Â£ coordinate line such that the spacing between each Â£ coordinate line is consistent with the prescribed grid cell volume. At the bounding rÂ¡ coordinate lines, the upstream axis of symmetry and the downstream boundary of Figure 5.2, the Â£ coordinate lines are required to intersect the coordinates at right angles. The projectile used in all the calculations is a 6-caliber, secant- ogive cylinder boattail configuration which is shown in Figure 5.4. There are experimentally determined pressure coefficients available at a series of points along the projectile surface for a range of Mach numbers in the transonic range [25]. The data will serve as a means of canparing the computed solutions. A grid network for this projectile problem with sting was generated using this technique with 70 points in the Â£ coordinate direction (IM=70) and 60 points in the rÂ¡ coordinate direction (JM=60). As shown in Figure 5.5, most of the points in the rÂ¡ coordinate direction were clustered near the projectile surface to resolve the boundary layer. The spacing at the projectile surface in the r? coordinate direction is 0.00002. The outer boundary is extended to four times the projectile length. In the streamwise direction, the points were clustered near the secant-ogive/cylinder and cylinder/boattail junctures in anticipation of the expansion waves that will occur in transonic flow. The first solution on this grid network is obtained with the Thin Layer code for Mach 0.96. The Reynolds number for this and all examples is 76,000. For all the results obtained, the time step used for a converged solution is 0.1. Figures 5.6 and 5.7 shew the convergence process of the solution algorithm for the indicated time spans. It is quite obvious that the solution oscillates in the region containing the 54 Figure 5.5 Initial fixed grid network for projectile with sting (70 x 60) 55 MACH 0.96 cPt i Figure 5.6 Converging solution for the Beam and Warming scheme Figure 5.7 Converged solution for the Beam and Warming scheme 56 Figure 5.8 Converged solution for the TVD scheme Figure 5.9 Least squared residue 57 shocks. However, after approximately 5000 time steps the solution does converge and agrees well with the experimental data as indicated by the comparison of the experimental and computed pressure coefficient along the projectile surface. Figure 5.8 shows the results obtained using the TVD code. The solution also agrees well with the experimental data; moreover it converged within 1600 time steps. This difference in the convergence ofthe two schemes is shown in Figure 5.9 which plots the least square residue R, versus the number of time steps N. Here the slower, oscillatory convergence of the Beam and Warming scheme is evident while the TVD scheme appears to converge monotonically. The projectile base flow problem is also solved on a fixed grid using the TVD code. The grid network for this calculation was obtained using a modified version of the technique of Steger et al. [16] in which the boundary condition on the downstream bounding rÂ¡ coordinate line has been changed so that the Â£ coordinate lines intersect with the x axis at right angles. The grid network obtained with IM=70 and JM=60 is shown in Figure 5.10. As can be seen, the Â£ coordinate lines now bend around the sharp comer, and end at the downstream axis of symmetry forming an O-type grid network. The computed pressure coefficient obtained for this problem using the TVD code is shown in Figure 5.11. Again, the solution converged in approximately 1600 time steps. It compares well, as before, over the secant-ogive and cylinder portions but then differs from the experimental results at the end of the boattail. No experimental data is available for the base region. There are a number of reasons that may cause the difference between the computed and experimental results near the comer. One reason may be the presence of a narrow sting used to support the experimental model that is not present in the computational 58 configuration. The difference could also be due to poor resolution in the projectile comer region or an inadequate turbulence model near the flow separation region. 59 Figure 5.10 Initial fixed grid network for projectile base flow problem (70 x 60) 60 Figure 5.11 Converged solution for projectile base flow problem using the TVD scheme CHAPTER VI PRELIMINARY TESTS AND MODIFICATIONS VI.1 The General Procedure Prior to coupling the Adaptive Grid code with the flow solvers to solve the projectile flow problem, it was found beneficial to modify the adaptive grid generation method. The following sections include an example showing the control over orthogonality and then describes two modifications made to the adaptive grid scheme to improve the general technique. It should be noted that each of the modifications is made to improve the method for the present flow problem and geometry and may not be necessary or useful in other problems. For each example the initial grid network used was that of Figure 5.5, or that of 5.10 if the base flow configuration is used in the example. In the adaptive grid code, the control function P was set equal to a constant so that the clustered points along each | coordinate line in the initial grid would tend to spread and equal spacing would result. The control function Q was defined using eqs. (4.3.2) with the distribution of eq. (4.3.3) so that the same point distribution along the r? coordinates in the initial grid would be retained in the grid obtained using the Adaptive Grid code. Each of the figures used in the examples show only the points in the reduced grid network. 61 62 VI.2 Orthogonality In order to demonstrate the effect of enforcing orthogonality, two cases were considered, one with the parameter A = 0.0, the second with A = 0.5. For the present projectile, the most difficult area to obtain an orthogonal grid is over the projectile nose. This is due to the sharp nose configuration which results in the two boundaries, the projectile surface and the upstream axis of symmetry, meeting at an obtuse angle. The two sets of coordinates, which run parallel to the boundaries, then tend to intersect at similar angles. This result is shown in Figure 6.1a, which shows an enlarged view of the grid in the nose section that resulted when A = 0.0. To improve the orthogonality in this region, another case for A = 0.5 was run, and as the results of Figure 6.1b indicate, the orthogonality of the grid coordinates near the nose was in general increased. For all the solutions to the transonic flow problem using the adaptive grid technique reported in Chapter VIII, A was set equal to 0.5. VI.3 Curvature One inherent result of grids obtained as solutions to equations based on the Laplacian operators is the effect of curved boundaries on the spacing of interior coordinate lines. In general, boundaries that are convex will attract coordinate lines and those that are concave will repel coordinate lines. In the context of the adaptive grid generation scheme presented here, this effect will alter the grid spacing dictated by the control functions, either increasing or decreasing the spacing that would result in the absence of curvature effects. The severity of the boundary curvature's influence on the grid is best seen in solving 63 Figure 6.1a Grid network using X=0.0 Figure 6.1b Grid network using \ =0.5 64 the adaptive grid generation equations for the base flow grid network. Recall that thel control function Q is defined so that the original grid point spacing in the normal direction would be retained. However, the point spacing resulting as a solution to the adaptive grid equations differs along each rÂ¡ coordinate line as shown in Figure 6.2a. The attraction of points is most prominent near the sharp comer at the projectile base where the effect is large due to the high curvature at the comer. A simple and effective technique to eliminate the effect of curvature can be developed by considering the solution of Laplace's equations between two concentric circles and C2 with radii R]_ and R2 respectively. V2Â£ = 0 (6.2.1) VrÂ¡ = 0 Solving Laplace's equation is equivalent to solving the adaptive grid equations with constant control functions and neglecting orthogonality (i.e. A = 0.0). In using constant control functions, we expect the coordinates to be equally spaced in each direction. It is advantageous to write the equations in polar coordinates (r,0), and it is sufficient here to seek a solution in which q is a function of 0 only and rÂ¡ is a function of r only. The Â£ coordinate lines then align themselves with the 0 coordinate direction and the rj coordinate lines became aligned with the radial coordinate r. The problem can then be formulated as 65 mw ui/i Figure 6.2a Grid network obtained without corrections for curvature Figure 6.2b Grid network obtained with corrections for curvature 66 d2Â£ = O de2 d2r? 1 dr? + =0 dr2 r dr (6.2.2) Â£ = 6 on Ci and C2 rÂ¡ = R]_ on C-]_ 1 = R2 on C2 The boundary conditions corresponds to uniform spacing along the Â£ coordinate lines coincident with the bounding concentric circles. The spacing along each coordinate direction is the inverse of the first derivatives and rÂ¡r which by integrating eqs. (6.2.2) are found to be to ~ a *?r (6.2.3) where a and b are constants of integration. These results indicate that the spacing in the Â£ direction is uniform while the spacing in the rÂ¡ direction is smaller near the inner circle and larger near the outer circle C2 The term in the governing equation responsible for this result, designated here as K, is K 1 d^ r dr (6.2.4) which can be viewed as the curvature of the intersecting Â£ coordinate line, (1/r), multiplied by the inverse of the spacing along the rÂ¡ 67 coordinate line. The effects of the curved f coordinates on the spacing along the r? coordinate lines can be eliminated for this problem simply by subtracting K from the governing equation for rÂ¡, v2Â£ = 0 (6.2.5) 0 1 dl Mzr, = 0 r dr This technique can be extended to the general adaptive grid equations by developing the curvature term K for each curvilinear coordinate line in the computational domain and then subtracting them from eqs. (2.2.7). The term needed to cancel the effects of a curved rÂ¡ coordinate line on the spacing along a f coordinate line is comprised of the curvature of the t) coordinate line and the spacing along the Â£ coordinate line. The curvature p of a rÂ¡ coordinate line is defined as P = | 7 | (6.2.6) ds where r is a unit tangent to the coordinate and is in the computational plane (xr? rYr]) Joe (6.2.7) where a is defined in eq. (2.3.5). By using the two relationships ds Sri dr? the curvature p can be shown to be 68 P =( Y V*r,r, a 3/2 (6.2.9) By multiplying this result with the inverse of the spacing along a Â£ coordinate (l//y), the following expression is obtained Ki = Yr)^-r)r] ~ ^-r)Yr]r] Jy 3/2 (6.2.10) which is the term that is subtracted from the adaptive grid equation for Â£ in order to cancel the effects of curved rÂ¡ coordinate lines. An expression for the influence of curved Â£ coordinates on the spacing along rÂ¡ coordinate lines can be obtained in a similar manner and is K2 = y^x^ xcyee Ja 7 3/2 (6.2.11) Thompson et al. [15] have derived similar expression, but do not implement them in the same way. It should be noted that in the computational plane, the terms and K2 contain second derivatives of the variables x and y, and when eqs. (2.2.7) are transformed onto the computational plane, the second order central difference approximations to these terms depend explicitly on the point (xifj, yi,j) However, in implementing the Newton-Raphson iterative solution procedure, the terms Ki and K2 are considered part of the source terms. Thus, the right hand side of eqs. (2.3.7) becomes 69 aJ P Si = ^ KiJ3 (6.2.12) s2 = qJ Q n Q k2J3 This is equivalent to neglecting the dependence of the finite difference approximation on the point (xj_f j, yifj) in obtaining eqs. (3.1.5) for the Newton-Raphson iterative procedure. The J3 term in eqs. (6.2.12) appears in the transformation of eqs. (2.2.7) into the computational plane. The adaptive grid generation equations, modified in this manner, will produce the grid shown in Figure 6.2b. In contrast to Figure 6.2a, it is clear that the effects of the curved boundaries have been successfully eliminated. VI.4 Internal Grid Boundaries One other problem that arose in the preliminary testing of the adaptive grid generation technique is the rapid upstream bending of the rÂ¡ coordinate lines emanating from the secant portion of the projectile surface. Figure 6.3a shows a grid network obtained using the adaptive grid generation equations with the curvature terms included and with A = 0.5. As can be seen, the rÂ¡ coordinate lines emanating from the secant portion of the projectile quickly bend forward to fill the upstream region of the domain which results in a very skew grid. This characteristic of the grid does not appear as a problem until a solution using the Beam and Warming scheme is sought on such a grid. Figure 6.3b shows a time history of the solution obtained on such a grid for Mach 0.96 and as can be seen the pressure coefficient oscillates in time over 70 Figure 6.3a Grid network upstream of the projectile Figure 6.3b Lack of convergence over secant ogive with Beam and Warming scheme 71 the secant portion of the projectile. The solution did not converge and instead the oscillations continued and eventually effected the solution downstream. In order to prevent the upstream bending of the r? coordinates, one rÂ¡ coordinate line in the initial grid network is designated as an internal boundary. This rÂ¡ coordinate line remains fixed in space during the iterativesolution procedure, thus no other rÂ¡ coordinate lines will pass through it. However, the points defining the internal boundary coordinate line are allowed to move along the coordinate so that they remain consistent with the spacing of the points along the adjacent rÂ¡ coordinate lines. This is accomplished by updating the arclength distribution along the internal boundary after each iterative sweep in which the new distribution is defined as the average of the arclength distributions along the two adjacent rÂ¡ coordinate lines. Figure 6.3c shews the grid obtained when the 12th rÂ¡ coordinate line from the projectile nose is designated as an internal boundary and as can be seen it prevents the rÂ¡ coordinate lines from bending upstream. It has been chosen to lie in a region in which no major grid clustering is expected so that it will not effect the adaptation of the grid network. The length and number of points used along the Â£ coordinate lines in each section may result in an abrupt change in the point spacing along the Â£ coordinate lines passing through the boundary as can be seen in Figure 6.3c. Therefore the control function P, controlling the spacing along the Â£ coordinate lines is locally modified near the internal boundary so that the grid point spacing varied smoothly through the internal boundary. This can be accomplished by requiring the normalized rate of change of grid point spacing along a Â£ coordinate line (s^/s^) 72 Figure 6.3c Grid network with internal boundaries Figure 6.4 Grid points near internal boundary 73 be identical at the two points adjacent to the internal boundary. Referring to Figure 6.4, if the h coordinate line designated as an internal boundary contains the i^*1 point along a Â£ coordinate line, we then require ^ . se . i-1 . se . i+1 (6.3.1) To implement this requirement into the solution procedure, consider the one dimensional equation for adaptation along a Â£ coordinate, eq. (2.3.8) which is written here as P fii (6.3.2) The left hand side of eq. (6.3.2) appears in the source terms of the two dimensional equations for adaptation, eqs. (2.3.7). Substitution of eq. (6.3.2)into the source terms for a grid point just ahead of the internal boundary results in the modified source terms at the point i-1 (si)i-l,j ai-l,jJi-l,j The modified source terms for the point just behind the internal boundary are obtained by replacing i-1 with i+1 in eq. (6.3.3). In order to apply the iterative solution technique a value for the terms (s^/s^) in eqs. (6.3.3)for each point adjacent to the internal boundary will be needed. Their final values are not known a priori since they will depend on the final position of the points in each section, therefore they are set to (6.3.3) i-l,i 74 the current value s e ii/j (si-2 ~2sj + si+2) (si+2-si-l)/2 (6.3.4) The grid that is obtained when using this technique is shown in Figure 6.5 in which it can be seen that the abrupt change in spacing across the internal boundary of the grid in Figure 5.3c has been reduced. The lose of internal boundaries provides a simple and effective way to control the grid points in the domain. In fact, it was found helpful to also designate the rÂ¡ coordinate line which separates the region over the sting from that over the projectile as an internal boundary to keep the rÂ¡ coordinate lines that originally began on the projectile surface from moving over the sting. 75 Figure 6.5 Smoothed grid spacing near the internal boundary CHAPTER VII THE SELF-ADAPTIVE CDMHJTATTONAL METHOD The adaptive grid generation technique can be incorporated naturally into the solution process of the flow solvers since the flow solver algorithms are time-marching. After a specified number of time steps, the grid network can be adapted to the control functions based on the current values of the flew variables. Thus the time-marching algorithm of the two flow solvers used here, the Thin Layer and TVD codes, proceeds as usual with intermittent calls to the Adaptive Grid code to update the grid network. In this approach, however, the relocation of the grid points in the physical domain due to the grid adaptation must be considered. The values of the flew variables at any time step are defined at the current grid point locations and any movement of the points will therefore distort the solution. This effect is critical for unsteady flow problems since the grid point movement will effect the solution accuracy. For steady flow problems, the grid point movement may effect the convergence rate and if it is large enough, it may cause the solution to diverge. There are two basic approaches to account for the grid point movement resulting when adapting the grid network. The first is to view the grid network as a moving grid and include the time metrics in the coordinate transformation of the Navier-Stokes equations. These metrics can be written in the computational plane as 76 77 t = Â£xxr + iy Yr 7t = 7xxr + 7yYr (7.1) The quantities x,. and yT can be approximated by using a backward difference expression xr (Xn+1-xn)/At Yr (yn+1-yn)/At two in which xn and yn are the grid point locations of the previous grid network and x1^1 are y1^1 are the values of the adapted grid network. This approach requires the grid to be adapted every time step which is inefficient for steady flow problems. Furthermore, it has been shown that special procedures must be used in evaluating the time derivative of the Jacobian [20] in order to maintain the conservative property of the finite difference schemes. Another approach, used here, is to interpolate the flow variables from the previous grid network onto the adapted grid network. This approach can be applied to both steady and unsteady flow problems and does not require the grid to be adapted every time step. Thus this approach can be used for both steady and unsteady flow problems. For the two- dimensional grid networks considered here, a linear interpolation is currently used. Each grid cell in the previous grid network is divided into two triangles. Once a grid point in the adapted grid network is located within a triangle of the previous grid network, the value of the flow variables at that point can be determined by interpolation. Referring to Figure 7.1, the value of a flow variable q at the point in the adapted grid network can be determined from the values of the flow 78 O previous grid network * current grid network Figure 7.1 Flow variable interpolation 79 variables at the vertices of the triangle as (qiAf + <%2a2 + <33 A3) q = (7.3) (Ai + A2 + A3) where A3, A2 and A3 are the areas of the triangles formed with the grid point in the adapted grid network. The basic algorithm for incorporating the adaptive grid generation technique into both the Thin Layer code and the TVD code consists of the following steps: 1. Input an initial grid network into the flow solver code and begin the time-marching algorithm. 2. Advance the solution a specified number of time steps, and submit the current grid network to the Adaptive Grid code. 3. Obtain an adapted grid network using the Adaptive Grid code. 4. Interpolate the flow variables onto the adapted grid network which new becomes the current grid network. 5. Repeat steps 2 through 5 until the solution converges for steady flew problems or repeat the steps until a specified time is reached for unsteady flow problems. This process which couples the flew solver and the adaptive grid generation technique is referred to here as the self-adaptive computational technique. For the two problems considered here, the axisymmetric transonic flew over a projectile with sting and the transonic projectile base flow problem, only steady flew solutions are sought. Since the adapted grid network depends on the flew variables via the control functions, the convergence of the solution implies also that the grid network ceases to change with time. In all the examples of the following chapters the grid networks of Figures 5.4 and 5.10 are used as 80 the initial grid. For the purpose of adapting the grid network, the source terras corrected to eliminate curvature effects, eqs. (6.2.12) were used in the adaptive grid generation equations. In solving the projectile problem with sting two internal r? coordinates were designated as internal boundaries. The first is over the secant ogive as shown in Figure 6.5, the second separates the grid network over the projectile from that over the sting and is located at approximately x=7.0. In solving the projectile base flow problem, no internal boundaries were used. CHAPTER VIII RESULTS AND DISCUSSION VIII.1 Choices for the Control Functions In applying the self-adaptive computational technique it is necessary to define the control functions for adaptation. For the transonic projectile problem considered here both shocks and expansion waves will occur approximately normal to the projectile surface and require refined grid point spacing in the Â£ coordinate direction for adequate resolution. The original choice for the function f in the control function P to obtain this clustering was the pressure gradient. However, it was found in the numerical experiments that the expansion waves occurring at the secant ogive/cylinder and cylinder/boattail junctures required smaller spacing than the shocks for good resolution. The pressure gradient was largest in the shocks and thus smaller spacing occurred in the shock regions rather at the expansion waves. A better choice for the control function was found to be the curvature of the pressure distribution along the streamwise coordinate direction, f = & 2 ( 1 + (** )2)3/2 as (8.1.1) 81 82 where the derivatives are with respect to the arclength s along each Â£ coordinate line. This function is largest at the expansion waves where the pressure gradient changed sign. At the base and peak of a shock eq. (8.1.1) also increased but it is not as large as the value near the expansion waves. Thus the smallest spacing will occur near the expansions and yet the grid points will also cluster, to a lesser extent, in the vicinity of shocks. In the normal direction, the important area requiring increased resolution is the boundary layer and the logical choice for the function g in the control function Q is the velocity gradient. However, it was found that the control function based on the velocity gradient resulted in too many points being placed in the boundary layer and led thus to a rapid stretching of the grid point spacing along the ry coordinate lines. Another choice for the function g, the exponential of the velocity gradient g = exp dn d s (8.1.2) where u is the velocity magnitude and s measures along the rÂ¡ coordinate line, was found to yield better results. Figure 8.1 shows a plot of the control function Q along a typical Sv coordinate line. The control function resulting from the exponential of the velocity gradient is similar to the control function corresponding to the grid point distribution in the fixed grid networks which is known to give good results provided enough points are used. The control function obtained using the velocity gradient however is large for many grid points and, 83 Q do6) Figure 8.1 Comparison of control functions Figure 8.2 Lack of convergence using a self-adaptive computational technique with Beam and Warming scheme 84 thus, results in too many points being clustered into the boundary layer region at the expense of adequate resolution in other areas. The minimum and maximum grid point spacing along the coordinate lines is 0.00002 and 6.0 respectively. In this direction it is important to have the small spacing at each point along the projectile surface in order to adequately resolve the viscous sublayer region. The parameter a of eq. (4.1.4) was determined for each rÂ¡ coordinate line so that the small spacing would occur for each point on the projectile surface. However, in the Â£ coordinate direction, the shocks and expansion waves are not expected to reach the outer boundary and the small spacing needed along the projectile surface is not required near that boundary. With the minimum and maximum grid point spacing prescribed as 0.02 and 0.22, respectively, a is determined only for the Â£ coordinate line coincident with the projectile surface and that value is used for all the remaining Â£ coordinate lines. This had the effect of increasing the minimum spacing as the magnitude of the shocks and expansion waves diminished near the outer boundary. VIII.2 Results Using the Self-Adaptive Computational Technique: The Beam and Warming Scheme The first application of the self-adaptive computational technique is with the Beam and Warming scheme for the calculation of the projectile flow with sting. The Mach number is 0.96 and the Reynolds number is 76,000, a case for which experimental data is available. The time step used in the calculation is 0.1 and the explicit and implicit dissipation terms are 2 and 4 times the time step respectively. The initial grid network used is that of Figure 5.5 with IM=70 and JM=60. In the first 85 attempts to use the adaptive grid technique, the grid network was adapted every 200 time steps. A series of calculated pressure coefficients obtained with this approach are shown in Figure 8.2. As can be seen, the pressure coefficient oscillates widely over the projectile surface and shows no indication of converging. In a series of attempts to obtain better results, the time step was decreased to 0.05, the dissipation terms were doubled and the number of time steps between adaptation were both increased to 500 and reduced to 50. However, similar results were obtained in each case. In the next case tried, the number of points in the streamwise direction was increased to 90 points, thus IM=90, and the grid network was adapted every time step. The results for this case are much better, however a converged solution was not obtained. Figure 8.3a shows the computed pressure coefficient at the times indicated, and as can be seen, the solution appears to have converged everywhere except in the shock boundary layer interaction regions. The shocks propagate upstream and downstream in a cyclic fashion. Figure 8.3b shows the same phenomenon occurring at a later time. In fact, the solution was continued for 8000 time steps and the propagation of the two shocks showed no indication of diminishing. Figure 8.4a shows the reduced grid network which resulted at 5000 time steps. As is evident, the grid point spacing in the rÂ¡ coordinate direction is clustered near the projectile surface with a distribution similar to that of the hyperbolic grid in Figure 5.4. The spacing varies somewhat in the shock and expansion wave regions, most likely as a result of enforcing orthogonality. The clustering in the streamwise direction 86 Figure 8.3a Computed solution between 4000 and 5000 time steps with Beam and Warming scheme Figure 8.3b Computed solution between 5000 and 6000 time steps with Beam and Warming scheme 87 Figure 8.4a Adapted grid network at 5000 time steps (90 x 60) Figure 8.4b Adapted grid network at 6000 time steps (90 x 60) 88 clearly indicates a response to the control function P, which is based on the curvature of the pressure function. Clustering occurs at both junctures on the projectile surface and the adaptation to the shock structures on the cylinder and boattail are also quite evident. Note that the smallest spacing appears to be located in the shocks well above the projectile surface. This occurs because there is more competition for the small grid spacing on the projectile surface since both the shocks and expansion waves are strong in that region. Figure 8.4b shows the grid network corresponding to a later time. In this Figure and those following, the Â£ coordinate lines are not shown in order to gain a better view of the grid spacing near the projectile surface. The grid configuration is basically the same as that of Figure 8.4a, however, the position of the shocks as indicated by the adaptation has changed. Ihe shock on the cylinder has propagated downstream and the boattail shock has moved upstream. Referring to the shock positions evident in Figures 8.3a and 8.3b it can be seen the clustering in the grid network is successfully following the shock patterns. In an attempt to improve these results the time step was reduced to 0.05 and the dissipation terms were again increased, however, similar results were obtained. Also the convergence criterion 5 of eq. (3.1.12) used in the iterative solution procedure for the grid generation equation was varied. At the value of 5=0.005, only one iteration was required for the grid network to converge. The rapid convergence occurred because the grid was updated every time step and the solution, and consequently the control functions, did not change much between grid adaptations. When the criterion was reduced to 6=0.0025, 1 to 3 iterations were required for 89 the grid to converge at each time step, however it had no effect on the cyclic propagation of the shocks. Although a converged solution was not obtained using the adaptive grid generation technique with the Beam and Warming scheme, there are still some encouraging results with respect to the adaptive grid generation procedure. The maximum and minimum grid point spacing occurring along the projectile surface in the adapted grid network were within 5 percent of their prescribed values. The average minimum spacing along all the r? coordinate lines was within 3 percent of the prescribed value of 0.00002 with a maximum deviation of 10 percent. The adapted grid networks shown in Figures 8.4a and 8.4b show that the adapted grid generation equations are capable of responding reliably and predictably to the chosen control functions and yield the prescribed maximum and minimum grid point spacing. The adaptive grid generation equations, thus, have been shown to provide a well adapted grid network for the transonic flow problems provided that a good choice is made for the control functions. However, the convergence of the solution was effected by the use of the adaptive grid generation procedure with the Beam and Warming scheme. One reason for the failure of the solution to converge may be due to interpolation. Error introduced due to interpolating the flow variables onto each successive adapted grid network may continuously perturb the solution, causing it to oscillate. Another cause of this phenomenon may be the form of the dissipation terms in the Beam and Warming scheme. The added dissipation terms contain derivatives with respect to the computational coordinates and were added to the scheme after the governing equations were transformed. Since these terms are a function of the transformed 90 coordinates any change in the grid network used will change the dissipation terms and consequently alter slightly the transformed governing equations. Thus the continuous adaptation of the grid network continuously alters the dissipation terms and may be the reason the solution does not converge. VIII.3 Results Using the Self-Adaptive Computational Technique: The TVD Scheme Since the solutions obtained using the Beam and Warming scheme did not converge, another scheme, the TVD scheme was also used. Recall that the primary difference between the two schemes is the dissipation terms added to the governing equations. The initial grid network used is that of Figure 5.4 with IM=70 and JM=60. The time step used is 0.1 and the grid was adapted every 200 time steps. For Mach 0.96, the solution converged easily within 2000 time steps. Only one iteration was needed for the grid network to converge at 1800 and 2000 time steps which also indicates a converged solution. Figure 8.5a shows the computed pressure coefficient using the adapted grid generation procedure compared to the solution obtained on the fixed grid network of Figure 5.4. The solutions agree well with experimental data over most of the projectile surface, except for the predicted position of the shock over the cylinder. The position predicted using the adaptive grid generation technique is centered on the cylinder whereas the position predicted using the fixed grid network is slightly upstream of the center. Figure 8.5b shows a shadowgraph for corresponding flow conditions in which it can be seen that the experimentally predicted shock location is centered over the cylinder portion of the projectile. Thus, the use of the adaptive grid generation technique improved the solution accuracy. Figure 8.5c shows 91 the adapted grid network corresponding to the converged solution. Again, clustering near the secant ogive/cylinder and cy1inder/boattail junctures is evident as well as the adaptation to the shocks. In comparison to the fixed grid network of Figure 5.4 it is evident that the grid spacing along the projectile surface is refined at the center of the cylinder portion due to the adaptation and is responsible for the increased solution accuracy. Two other cases were run using the identical procedure used for the case of Mach 0.96. Figure 8.6a shows a comparison of the solution obtained using both the fixed grid network and the adaptive grid generation technique for Mach 0.91, another case for which experimental data is available. The only improvement obtained lasing the adaptive grid generation procedure is an increase in the peak of the shock occurring over the cylinder portion of the projectile. However, the peak still does not reach the experimentally predicted value. In an attempt to increase the accuracy in that region, 10 more points were added in the Â£ coordinate direction, thus IM=80, however no further improvement in the solution obtained with the adaptive grid generation technique was evident. The results thus may be the best that the TVD scheme can yield given the current grid configuration and the shortened pressure peak is probably due to the larger dissipative error associated with the TVD scheme. Figure 8.6b shows the grid network corresponding to the converged solution. Again grid point clustering in response to the curvature of the pressure is evident, and as can be seen by the adaptation of the grid network, both shocks occur directly downstream of the expansion waves. 92 Figure 8.5a Comparison of computed solutions on fixed and adaptive grid networks with TVD scheme F Figure 8.5b SOCBT shadowgraph for Mach 0.96 (reprinted from reference 24) 93 Figure 8.5c Adapted grid network for converged solution for Mach 0.96 Figure 8.6a Comparison of computed solutions on fixed and adapted grid networks with TVD scheme 94 Figure 8.6b Adapted grid network for converged solution for Mach 0.91 Figure 8.7a Comparison of computed solutions on fixed and adapted grid networks with TVD scheme 95 Figure 8.7b Adapted grid network for converged solution for Mach 1.2 Figure 8.8 Setting maximum value for the control function 96 The results for a third case, Mach 1.2, are shewn in Figures 8.7a and 8.7b. There is no noticeable difference between the computed results on the fixed grid network and the results obtained using the adaptive grid generation procedure except near the projectile nose where the grid configuration differs the most. In the adapted grid network, clustering again is evident near the junctures, however no shocks are present for this case, a result also evident in the grid network. The equal accuracy of the solutions obtained with the adapted grid generation procedure and on the fixed grid network is a result of the similarity of the grid networks. In comparing the adapted grid networks of Figures 8.6b and 8.7b to the fixed grid network of Figure 5.4 it is evident that they have basically the same features in that the grid point clustering along the Â£ coordinates is concentrated near the two junctures. Thus the original fixed grid network is well adapted to the problem since the shocks for case of Mach 0.91 are so close to the expansion waves and no shocks occur in the case of Mach 1.2. However, it should be noted that the adaptive grid generation procedure successfully produced a good adapted grid network without any knowledge of the position or orientation of the important flew phenomenon. The choice and construction of the proper control functions of course are important in obtaining these results. In the context of these results it can be concluded that the adaptive grid generation procedure can reliably produce a good adapted grid network provided good choices are made for the control functions. 97 VIII.4 Application to the Projectile Base Flow Problem In a final case, the self-adaptive computational technique was applied to the projectile base flow problem. The flow conditions were Mach 0.96 and a Reynolds number of 76,000. As in the applications to the projectile with sting, the time step used in the TVD scheme is 0.1 and the grid is adapted every 200 time steps. The grid network of Figure 5.10 was used as the initial grid, however in the adaptive grid procedure no internal boundaries were used. The coordinate lines emanating from the boattail are expected to bend downstream in the current configuration and, thus, should prevent the lines emanating from the secant ogive from bending to far foward. Before applying the self-adaptive computational technique to compute the solution, the 70x60 grid network was adapted to the converged solution obtained on the 70x60 fixed grid network. It was found that the control function in a small region near the base comer was so large that most of the grid points clustered in that region and resulted in poor resolution near the shocks and expansion waves. The method used to alleviate this problem, is to prescribe a maximum allowable value for the curvature of the pressure distribution. Any value of the curvature above this value was simply set to the prescribed maximum value. This approach is depicted in Figure 8.8. It introduces another parameter but results in a much better grid network and makes the self-adaptive grid generation much more reliable. Currently the maximum value used is 30, which was determined after some experimentation. The calculated pressure coefficent obtained using the self-adaptive computational technique is shown in Figure 8.9a; the distribution over 98 the base region is shown in Figure 8.9b. There is generally good agreement between the two calculated solutions with the experimental data over the projectile surface with a slight discrepency near the shock on the cylinder. At the base comer and along the base the two calculated solutions are different though their trends are similar. Without any experimental data for comparison it is difficult to judge the accuracy. Furthermore, the thin layer approximation and the current turbulence model may not be valid near the base since a separated wake region exists. The adaptive grid network corresponding to the converged solution is shown in Figure 8.10a and an enlarged view of the grid network near the base comer is shown in Figure 8.10b. In comparing this grid network to that of Figure 8.5c for the projectile with sting problem it is evident that the r? coordinate lines emanating from the front portion of the projectile do bend upstream, but not so drastically, and do not prevent the convergence of the solution using the TVD scheme. There is the same general adaptation to the expansion waves and shocks, however the adaptation to the shock over the boattail appears to bend downstream, a result due to the smoothing of the control functions in the computational plane. Adaptation of the grid network near the base comer region is also evident. Note that the grid lines remained nearly orthogonal in the base region even though the normal direction changes direction rapidly near the sharp comer. In future work, the orthogonlity near the comer may be improved by decreasing the spacing of the Â£ coordinate lines in the reduced grid network. Also, the grid spacing between the rÂ¡ coordinate lines is fairly uniform and smooth and the coordinate lines did not 99 overlap at the comer, a result attributed to the elliptic nature of the adaptive grid generation equations. The adaptive grid generation procedure, thus, can successfully provide a good adapted grid for this complex geometry. 100 0. 0. -0. -0. 2 0 a 2 4 / a - MACH 0.96 b x s X\ f .-..V. - adapted grid fixed grid b i 0.37 0.0 Figure 8.9 Comparison of computed solutions for the projectile base flow problem with TVD scheme 101 Figure 8.10 Adapted grid network for the converged solution for Mach 0.96 CHAPTER IX CONCLUDING REMARKS The self-adaptive computational technique developed here has been shown to provide good adaptive grid networks for the transonic projectile flow problems considered. The equations governing the grid adaptation were developed, based on a variational approach, to provide adaptation independently in each coordinate direction and also enhance both grid orthogonality and smoothness. The resulting set of elliptic partial differential equations provide explicit control of orthogonality and adaptation while smoothness is inherent in the ellipticity of the governing equations. Smoothness can also be controlled in the definition of the control functions. A series of modifications were made to improve both the adaptive grid generation equations and their solution procedure. local scaling of the equations was found useful for highly stretched meshes, and an effective technique to eliminate the effects of curved boundaries, especially useful in the projectile base flow problem, was developed. Also, the numerical solution procedure for the adaptive grid generation equations was improved to increase the poor convergence rate due to the high aspect ratio grid cells present in the grid networks used for transonic viscous flow problems. Finally, the use of internal grid boundaries, which were used in the projectile flow problem with sting, helped to provide a better grid network by restricting the movement of the grid points. 102 103 In coupling the adaptive grid generation procedure with the flow solvers it was shown that the adaptive grid generation equations could provide good adapted grid networks for the transonic projectile flow problems. It was found that the control function based on the curvature of the pressure distribution was more effective than the pressure gradient when both shocks and expansion waves are present. Also, in using a control function based on the velocity gradient, it was shown that the adaptive grid generation procedure could provide the extremely refined grid point spacing necessary in the boundary layer region. The most important feature of the adaptive grid generation procedure developed here is the simple relationship between the control functions and the grid point spacing which is maintained in an elliptic set of partial differential equations. As more complex problems are attempted, the control functions for each coordinate direction can be carefully designed, based on various physical gradients, prescribed grid point distributions and combinations thereof, to obtain the desired adaptation. The elliptic nature of the adaptive grid generation equations will be helpful, especially for complex geometries, in maintaining a one to one mapping and preventing highly distorted grid networks from occuring. frOT %Aux uA$x = Â£ zp i* + ^ ^ x = Mi ^ 5aJx + i^i! |x + ueiS )xUxz ii| Sx = >tXq Ze 2e er/((A)xxT^x + (x)xxT A-) = xxk Er/((A)^xq5x + (x)^xt A-) = ^XU Er/((A)^^qx + (x)^^q A-) = xxu er/( (A)xxq^x + (x)xxq **A-) = Er/( (AJ^q^x + (x)^xq UA~) = ^x5 cr/((AJ^^I^X + (x)AAq ^A-) = xx AXrj XIQNaddV REFERENCES 1. Kutler, P, "A Perspective of Theoretical and Applied Computational Fluid Dynamics," AIAA Journal, Vol. 23, NO. 3, pp. 328-341, 1984. 2. Thompson, J.F., Thames, F.C., and Mastin, C.W., "Automated Numerical Generation of Body-Fitted Curvilinear Coordinate System for Field Containing Any Number of Arbitrary Two-Dimensional Bodies," J. Comp. Phys., Vol. 15, pp. 299-319, 1974. 3. Mastin, C.W., "Error Induced By Coordinate Systems," Numerical Grid Generation. Ed. J.F. Thompson, North-Holiand, p. 41, 1981. 4. Hsu, C.C., and Chang, T.H., "On a Numerical Solution of Incompressible Turbulent Boundary Layer Flow," Finite Element Flow Analysis. Ed. T. Kawai, Univ. of Tokyo Press, Tokyo, pp 219-226, 1982. 5. Saltzman, J. and Brackbill, J.U., "Application and Generalization of Variational Methods for Generating Adaptive Meshes," Numerical Grid Generation. Ed. J.F. Thompson, North-Holland, pp. 865-884, 1981. 6. Hsu, C.C. and Tu, C.G., "An Adaptive Grid Generation Technique Based on Variational Principals for Transonic Projectile Aerodynamic Calculation," To Appear in Int. J. Numerical Methods in Fluids. 7. Nakahashi, K., and Deiwert, G.S., "A Self-Adaptive Grid Method with Applications to Airfoil FLow," AIAA Paper 85-1525. 8. Dywer, H.A., Smooke, M.D. and Kee, R.J., "Adaptive Gridding for Finite Difference Solutions to Heat and Mass Transfer Problems," Numerical Grid Generation. Ed. J.F. Thompson, North-Holland, New York, p. 339, 1982. 9. Gnoffo, P. A., "A Vectorized, Finite Volume, Adaptive-Grid algorithm for Navier-Stokes Calculations," Numerical Grid Generation. Ed. J. F. Thompson, North-Holland, New York, p. 819, 1981. 10. Thompson, J.F., "A Survey of Dynamically-Adaptive Grids In the Numerical Solution of Partial Differential Equations," Applied Numerical Mathematics 1. North-Holland, New York, 1985. 11. Rai, M.M. and Anderson, D.A., "Grid Evolution in Time Asymptotic Problems," J. of Comp. Rhys., Vol 43, pp. 327-344, 1981. 105 106 12. Ames, W.F., "Numerical Methods for Partial Differential Equations.11 Academic Press, Inc., New York, 1977. 13. Tu, C.G., "Development of An Adaptive Grid Generation Code for Projectile Aerodynamic Computation," Ph.D. Dissertation, University of Florida, 1985. 14. Pierson, B.L. and Kutler, P., "Optimal Nodal Point Distribution for Improved Accuracy in Computational Fluid Dynamics," AIAA Journal, Vol. 18, NO. 1, pp. 49-54, 1980. 15. Thompson, J. F., Warsi, Z. U. A. and Mast in, C. W. Numerical Grid Generation: Foundations and Applications. North-Hoiland, New York, pp. 215-221, 1985. 16. Nietubicz, C.J., Heavey, K.R. and Steger, J. L., "Grid Generation Techniques for Projectile Configurations," Proc. 1982 Army Numerical Analysis and Computers Conference, ARO Rept. 82-3, February 1982. 17. Warming, R.F. and Beam, R., "On the Construction and Application of Implicit Factored Schemes for Conservation Laws," SIAM-AMS Proceedings, Vol. 2, pp. 85-99, 1978. 18. Pulliam, T.H., and Steger, J.L., "Implicit Finite-Difference Simulations of Three-Dimensional Compressible Flow," AIAA Journal, Vol. 18, NO. 2, pp. 167-169, 1980. 19. Nietubicz, C.J., Pulliam, T.H. and Steger, J.L., "Numerical Solutions of the Azimuthal Invariant Thin Layer Navier-Stokes Equations," AIAA Paper 79-0010, Presented at AIAA 17th Aerospace Sciences Meeting, New Orleans, Louisiana. 20. Anderson, D.A., Tannehill, J.C. and Pletcher, R.H., Computational Fluid Mechanics and Heat Transfer. McGraw-Hill, New York, p. 546, 1984. 21. Baldwin, B.S., and Lomax, H., "Thin Layer Approximation and Algebraic Model for Separated Turbulent Flows," AIAA Paper 78-257, Presented at AIAA 16th Aerospace Sciences Meeting, Huntsville, Alabama, 1978. 22. Yee, H.C., "Linearized Form of Implicit TVD Schemes for the Multidimensional Euler and Navier-Stokes Equations," Comp. & Math, with Appls., Vol. 12a, Nos. 4/5, pp. 413-412, 1986. 23. Shiau, N. and Hsu, C.C., "A Diagonalized TVD Scheme for Turbulent Transonic Projectile Aerodynamic Calulations," submitted to 8th Int. Conf. on Computing Methods in Applied Sciences and Engineering, Versailles, France. 24. Sahu, J., Nietubicz, C.J. and Steger, J.L., "Numerical Computation of Base Flow for a Projectile at Transonic Speeds," ARERL-TR-02459, U.S. Army Ballistic Research Laboratory, June 1983. 107 25. Kayser, L.D. and Whiton, F., "Surface Pressure Measurements on a Boattailed Projectile Shape at Transonic Speeds," ARBRL-MR-03161, U.S. Army Ballistic Research laboratory, March, 1982. BIOGRAPHICAL SKETCH Christopher William Reed was bom in Rahway, New Jersey, in 1960. He graduated from the Bolles School in Jacksonville, Florida, in 1978 and received a B.S. in engineering science and mechanics from the Georgia Institute of Technology. He earned his M.S. from the University of Florida Department of Engineering Sciences in 1984 and in 1986 was married to Shelley Anne Baylis. He received his Fh.D. from the University of Florida Department of Engineering Sciences in 1987. 108 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 Fhilosophy. Chen-Chi Hsu, Chairman Professor of Engineering Sciences 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 (fegree of Doctor of Philosophy. Bernard tadori Professor of Engineering Sciences 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. Z//tn cl ?/ ')<'<* Ulrich H. Kurzweg Professor of Engineering Sciences 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 Fhilosophy. David W. Mikolaitis Assistant Professor of Engineering Sciences 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. C* ll^bz David C. Wilson Associate 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. August 1987 liaAcjf &. Dean, College of Engineering Dean, Graduate School 53 coordinate direction, defining each successive Â£ coordinate line such that the spacing between each Â£ coordinate line is consistent with the prescribed grid cell volume. At the bounding rÂ¡ coordinate lines, the upstream axis of symmetry and the downstream boundary of Figure 5.2, the Â£ coordinate lines are required to intersect the coordinates at right angles. The projectile used in all the calculations is a 6-caliber, secant- ogive cylinder boattail configuration which is shown in Figure 5.4. There are experimentally determined pressure coefficients available at a series of points along the projectile surface for a range of Mach numbers in the transonic range [25]. The data will serve as a means of canparing the computed solutions. A grid network for this projectile problem with sting was generated using this technique with 70 points in the Â£ coordinate direction (IM=70) and 60 points in the rÂ¡ coordinate direction (JM=60). As shown in Figure 5.5, most of the points in the rÂ¡ coordinate direction were clustered near the projectile surface to resolve the boundary layer. The spacing at the projectile surface in the r? coordinate direction is 0.00002. The outer boundary is extended to four times the projectile length. In the streamwise direction, the points were clustered near the secant-ogive/cylinder and cylinder/boattail junctures in anticipation of the expansion waves that will occur in transonic flow. The first solution on this grid network is obtained with the Thin Layer code for Mach 0.96. The Reynolds number for this and all examples is 76,000. For all the results obtained, the time step used for a converged solution is 0.1. Figures 5.6 and 5.7 shew the convergence process of the solution algorithm for the indicated time spans. It is quite obvious that the solution oscillates in the region containing the 11 (2.1.3) into the governing equations makes Â£ and rÂ¡ the independent spatial variables and now a numerical algorithm can be developed using finite difference approximations on the equally spaced computational domain. Since the metric coefficients depend only on the mapping, virtually any domain may be mapped into the computational domain and solved using the numerical algorithm. The metric coefficients can be evaluated for a numerically generated transformation by considering the metrics of the inverse transformation. The relationship between the metric and inverse metrics is, for two dimensions, Â£ x ?y 1 ' Yr, Y ' >1 >1 * J . ~xn . where J is the Jacobian of the transformation, J = x^Yr, ~ xnYt (2.1.5) The derivation of these equations can be found in reference [13]. Since Â£ and rÂ¡ are the independent variables for the inverse metrics, the inverse metrics can be evaluated numerically in the computational domain and subsequently used to determine the metric coefficients. The primary disadvantage of using boundary fitted curvilinear coordinate systems is the appearance of additional truncation error terms in the finite difference approximations on the computational plane. Mastin [3] has shown, for instance, that the truncation error T for a second order finite difference approximation of the first derivative is, 54 Figure 5.5 Initial fixed grid network for projectile with sting (70 x 60) 42 2 i=2,IM-1 j=l,JM (4.2.1) 2 i) i=l,IM j=2,JM-1 where the value of w used here is 0.4. The control function at a point is replaced by a weighted average of itself and the two adjacent points along a coordinate line. Currently, two consecutive sweeps of this algorithm are made over the domain for each of the control functions, P and Q. In this algorithm there is no accounting for the varying spacing between points and thus this smoothing is in the computational plane. IV.3 Other Control Functions In many instances a certain point distribution in one coordinate direction is known to yield good results and it would be advantageous to maintain this prescribed distribution in that direction while allowing adaptation to the solution in the other directions. For instance, in some of the viscous transonic flow problems solved here, it is known that a distribution in the direction normal to the surface based on an exponential function will yield good results provided enough points are used. A procedure to incorporate specified point distributions into the adaptive grid scheme is easily developed by using eq. (2.3.10). Consider some given grid point distribution s in the rÂ¡ coordinate direction which is to be reproduced by the adaptive grid generation equations 52 Figure 5.3 Boundary configuration for projectile base flow problem SOCBT Figure 5.4 SOCBT projectile configuration 46 Here, p is the density, p is the pressure and E is the total specific energy. The velocities u,v and w correspond to the x, y and z directions, respectively, as indicated in Figure 5.1 and U, V, and W are the contravariant velocity components corresponding to the Â£, rÂ¡ and f directions respectively. The coefficient of viscosity is p., Pr is the Prandtl number, Re is the Reynolds number, 7 is the ratio of specific heats and a is the local speed of sound. The details of the transformation and application of the thin layer approximation to obtain eqs. (5.1.1) are available in reference [13]. The assumption of a calorically and thermally perfect gas was made in closing the system of equations. The coefficient of viscosity pi is considered variable and it is also necessary to account for the effects of turbulence for the present applications. Under the assumption of f invariance, only a two dimensional grid network containing Â£ and r? coordinates in a plane through the axis of symmetry will be required. It should be noted, however, that the metric coefficients represent the three dimensional coordinate transformation and are therefore different from the expressions of eqs. (2.1.4) and the Jacobian J* now represents the grid cell volume. V.2 Solution Algorithms for the Governing Equations The implicit spatially factored numerical scheme of Beam and Warming [17] is used to solve the transformed governing equations. The scheme is implemented for three dimensional problems in a code due to Steger and Pulliam [18]; the code used here is a modified version which was developed by Nietubicz et al. [19] for the efficient calculation of axisymmetric projectile problems and is referred to here as the Thin 29 Axi,j = 2P 2 ^r? 1 Ayi^ = V 2Q with that for a cell of large aspect ratio AR 1, mj.2 (-L-, L'-1 2P 1+AR (3.2.4) _1 AR . Ayi -i ( ) -L,J 2Q 1+AR (3.2.5) it becomes evident that the reduction in the grid point motion in the x direction for large AR is proportional to (1/AR). For aspect ratios of the order (105) this result becomes prohibitive, rendering the iterative procedure extremely inefficient. An example has been formulated to demonstrate this deficiency. In an initial grid network, as shewn in Figure 3.2, the points along the Â£ coordinates are clustered in one region and the points along each rÂ¡ coordinate are clustered near the bottom surface to form the highly stretched spacing typical of grid networks used in transonic flow calculations. The spacing at the surface is approximately 105 times that in the upper region of the grid. In the adaptive grid generation equations, the control function P is set equal to a constant so that the points will redistribute themselves with equal spacing along the Â£ coordinates. Using a technique described in Chapter IV, the control function Q is defined so that the highly stretched spacing along each rÂ¡ 28 coordinates. For this example the aspect ratio AR is equal to (1/b). Assume for purposes of demonstration that the two terms (Qv/Q) and (P^/P) are non-zero at the center point and will therefore cause the point to move. We can obtain the changes in the point location by evaluating eqs. (3.1.7) for this simple case. Many terms for this example are zero xr?=o yz=0 y?7=0 y^=0 x*e= y^=0 xr,r,=0 Yr,r,=0 The coefficients became a3=b3 b]_=0 c3=0 d3=b2 a2=0 b2=0 c2=0 d2=0 a3=b b3=0 c3=0 d3=l (3.2.1) (3.2.2) and upon substitution into eqs. (3.1.7) yield Axi, j ( JL 2P l+b AYi,j = 2Q l+b (3.2.3) Writing these in terms of the aspect ratio AR, its effect on the grid point motion becomes clear. If we compare the solution for a square cell , i.e. AR=1, 9 mesh in the computational domain. Furthermore, once a successful finite difference algorithm is developed on the computational domain, any physical domain may be mapped into this region and solved with the same numerical algorithm. The task of numerical grid generation is to define the transformation between the Cartesian and computational coordinate systems which can be viewed as the development of a curvilinear coordinate system in the physical domain. The general mapping between the two coordinate systems is = (x,y) (2.1.1) 1 = v (x,y) Owing to the discrete nature of finite difference approximations, it is sufficient to define a finite set of coordinate lines, the intersections of which define the grid points in the physical domain. Thus for each grid point (x,y) in the physical domain, there are two curvilinear coordinates (Â£,??) that intersect at that point and designate its location. It is convenient to view this mapping in the computational domain, where for each pair (|,r?) there corresponds a point (x,y) in the physical domain. This correspondence can be conveniently expressed as the inverse mapping x = x (Â£ rÂ¡) Y = (2.1.2) 47 Figure 5.1 Coordinate systems for projectile configuration B Figure 5.2 Boundary configuration for projectile with sting 67 coordinate line. The effects of the curved f coordinates on the spacing along the r? coordinate lines can be eliminated for this problem simply by subtracting K from the governing equation for rÂ¡, v2Â£ = 0 (6.2.5) 0 1 dl Mzr, = 0 r dr This technique can be extended to the general adaptive grid equations by developing the curvature term K for each curvilinear coordinate line in the computational domain and then subtracting them from eqs. (2.2.7). The term needed to cancel the effects of a curved rÂ¡ coordinate line on the spacing along a f coordinate line is comprised of the curvature of the t) coordinate line and the spacing along the Â£ coordinate line. The curvature p of a rÂ¡ coordinate line is defined as P = | 7 | (6.2.6) ds where r is a unit tangent to the coordinate and is in the computational plane (xr? rYr]) Joe (6.2.7) where a is defined in eq. (2.3.5). By using the two relationships ds Sri dr? the curvature p can be shown to be 97 VIII.4 Application to the Projectile Base Flow Problem In a final case, the self-adaptive computational technique was applied to the projectile base flow problem. The flow conditions were Mach 0.96 and a Reynolds number of 76,000. As in the applications to the projectile with sting, the time step used in the TVD scheme is 0.1 and the grid is adapted every 200 time steps. The grid network of Figure 5.10 was used as the initial grid, however in the adaptive grid procedure no internal boundaries were used. The coordinate lines emanating from the boattail are expected to bend downstream in the current configuration and, thus, should prevent the lines emanating from the secant ogive from bending to far foward. Before applying the self-adaptive computational technique to compute the solution, the 70x60 grid network was adapted to the converged solution obtained on the 70x60 fixed grid network. It was found that the control function in a small region near the base comer was so large that most of the grid points clustered in that region and resulted in poor resolution near the shocks and expansion waves. The method used to alleviate this problem, is to prescribe a maximum allowable value for the curvature of the pressure distribution. Any value of the curvature above this value was simply set to the prescribed maximum value. This approach is depicted in Figure 8.8. It introduces another parameter but results in a much better grid network and makes the self-adaptive grid generation much more reliable. Currently the maximum value used is 30, which was determined after some experimentation. The calculated pressure coefficent obtained using the self-adaptive computational technique is shown in Figure 8.9a; the distribution over 93 Figure 8.5c Adapted grid network for converged solution for Mach 0.96 Figure 8.6a Comparison of computed solutions on fixed and adapted grid networks with TVD scheme 4 resolved a separated shear layer by adapting to the velocity gradient. One dimensional adaptation can be extended to higher dimensions by successive adaptation in each coordinate direction. Nakahashi and Deiwert [7] have used the spring analogy in such an approach and added torsional springs that connect intersecting coordinate lines to control orthogonality. They solved a transonic viscous airfoil problem in which the grid was adapted to the density gradient in the streamwise direction to resolve shocks and adapted to the velocity gradient normal to the airfoil surface to resolve the boundary layer. The one-dimensional approach to multi-dimension adaptation has the advantage of efficiency since the grid can usually be obtained using a marching solution algorithm. Another advantage is the independence of adaptation in each coordinate direction. However, it is difficult to maintain smoothness in such approaches, and highly skewed grids may result [10]. Rai and Anderson [11] have used a different approach for two dimensional adaptation. In their scheme each grid point induces a velocity, either repulsive or attractive, on every other point. By choosing some flow gradient to indicate the need for point clustering, each point with a gradient larger than the average gradient attracts points, and those with values below the average repel points. By summing the induced velocities from each point, the velocity of each point is determined and can be integrated over a time step to determine the new point location. In applying this scheme to the two dimensional linearized viscous Burger's equation, they adapted the grid to the derivative of the dependent variable in each direction and showed that such an approach will increase solution accuracy. 78 O previous grid network * current grid network Figure 7.1 Flow variable interpolation CHAPTER IV CONTROL FUNCTIONS IV. 1 The Basic Form It is the role of the control functions in the equations governing grid adaptation to dictate the need for increased resolution and, as indicated in section II. 1, they should be judicially chosen so that the resulting grid resolution reduces what would otherwise be large truncation error. For instance, Pierson and Kutler [14] chose for the control function, when solving a one dimensional inviscid Burger's equation, the third derivative of the dependent variable which is also the leading truncation error term in their finite difference approximation to Burger's equation. Their one dimensional adaptive grid scheme was thus an attempt to minimize a measure of the truncation error over the domain. In more complex problems, however, with more dependent variables and in higher dimensions, the truncation error expressions are complex and cumbersome to calculate. Therefore, in practice, an understanding of the physical problem and experience guide the choices for the control functions. The general control function considered here is P = 1 + 7pf Q = 1 + 7qg (4.1.1) 39 17 equations that govern the grid coordinates in the physical domain. II.3 Inversion of the Grid Generation Equations In order to solve the grid generation equations, they are transformed onto the curvilinear coordinate system. This is done by inverting eqs. (2.2.7) and (2.2.9) to make x, y, and s the dependent variables. The curvilinear coordinates then become the independent variables and the grid generation equations are solved on the computational domain which is equivalent to defining the inverse mapping of eqs. (2.1.2). All the advantages of using the computational domain to solve the governing equations discussed previously will be realized here as well. There are two approaches to inverting the equations. One is to transform the functional Ip onto the computational domain and then apply the Euler-Lagrange equations for x and y. Another approach, used here is to derive the corresponding expressions on the computational domain for each term in eqs. (2.2.7) and (2.2.9) and then to substitute them into the equations. The first term in eq. (2.2.7) can be written using eq. (2.1.3) as Â£xx = (Â£x + 7X ~ ) Â£x (2.3.1) dr) Using eq. (2.1.4) to transform the metrics this expression becomes Yr, d Yz d Yv Â£ xx -( ) (2.3.2) J J dr) J and after differentiating, the relationship is found to be 74 the current value s e ii/j (si-2 ~2sj + si+2) (si+2-si-l)/2 (6.3.4) The grid that is obtained when using this technique is shown in Figure 6.5 in which it can be seen that the abrupt change in spacing across the internal boundary of the grid in Figure 5.3c has been reduced. The lose of internal boundaries provides a simple and effective way to control the grid points in the domain. In fact, it was found helpful to also designate the rÂ¡ coordinate line which separates the region over the sting from that over the projectile as an internal boundary to keep the rÂ¡ coordinate lines that originally began on the projectile surface from moving over the sting. 77 t = Â£xxr + iy Yr 7t = 7xxr + 7yYr (7.1) The quantities x,. and yT can be approximated by using a backward difference expression xr (Xn+1-xn)/At Yr (yn+1-yn)/At two in which xn and yn are the grid point locations of the previous grid network and x1^1 are y1^1 are the values of the adapted grid network. This approach requires the grid to be adapted every time step which is inefficient for steady flow problems. Furthermore, it has been shown that special procedures must be used in evaluating the time derivative of the Jacobian [20] in order to maintain the conservative property of the finite difference schemes. Another approach, used here, is to interpolate the flow variables from the previous grid network onto the adapted grid network. This approach can be applied to both steady and unsteady flow problems and does not require the grid to be adapted every time step. Thus this approach can be used for both steady and unsteady flow problems. For the two- dimensional grid networks considered here, a linear interpolation is currently used. Each grid cell in the previous grid network is divided into two triangles. Once a grid point in the adapted grid network is located within a triangle of the previous grid network, the value of the flow variables at that point can be determined by interpolation. Referring to Figure 7.1, the value of a flow variable q at the point in the adapted grid network can be determined from the values of the flow BIOGRAPHICAL SKETCH Christopher William Reed was bom in Rahway, New Jersey, in 1960. He graduated from the Bolles School in Jacksonville, Florida, in 1978 and received a B.S. in engineering science and mechanics from the Georgia Institute of Technology. He earned his M.S. from the University of Florida Department of Engineering Sciences in 1984 and in 1986 was married to Shelley Anne Baylis. He received his Fh.D. from the University of Florida Department of Engineering Sciences in 1987. 108 CHAPTER II THE ADAPTIVE GRID GENERATION EQUATIONS II. 1 The Curvilinear Coordinate System The use of boundary fitted curvilinear coordinate systems in solving equations on arbitrarily shaped domains was originally motivated due to problems in applying finite difference approximations at curved boundaries. Because curved boundaries do not in general coincide with the rectilinear Cartesian coordinates, the grid points defined along these coordinates do not usually lie on the bounding curve. It therefore becomes difficult to accurately represent the boundaries [2] and special procedures must be used in those regions [12]. A solution to this problem, which has become an effective and versatile approach, is to map the arbitrarily shaped physical domain (x,y) into a rectangular computational domain (Â£, r?). Such a mapping for a typical projectile flow domain is shown in Figure 2.1. As indicated, the curved Â£ and rÂ¡ coordinate lines coincident with the boundaries in the physical domain map onto the edges of the rectangular computational domain. The boundary grid points, being defined by the coordinate intersections, are coincident with the boundary in the computational plane and thus no special treatment for the boundary conditions is required. Another advantage is that virtually any set of nonuniformally spaced curvilinear coordinates in the physical domain map into an equally spaced regular 7 VIII RESULTS AND DISCUSSION 81 VIII.1 Choices for the Control Functions 81 VIII.2 Results Using the Self-Adaptive Computational Technique: The Beam and Warming Scheme 84 VIII.3 Results Using the Self-Adaptive Computational Technique: The TVD Scheme 90 VIII.4 Applications to the Projectile Base Flow Problem 97 IX OCNCIUDING REMARKS 102 APPENDIX 104 REFERENCES 105 BIOGRAPHICAL SKETCH 108 V 91 the adapted grid network corresponding to the converged solution. Again, clustering near the secant ogive/cylinder and cy1inder/boattail junctures is evident as well as the adaptation to the shocks. In comparison to the fixed grid network of Figure 5.4 it is evident that the grid spacing along the projectile surface is refined at the center of the cylinder portion due to the adaptation and is responsible for the increased solution accuracy. Two other cases were run using the identical procedure used for the case of Mach 0.96. Figure 8.6a shows a comparison of the solution obtained using both the fixed grid network and the adaptive grid generation technique for Mach 0.91, another case for which experimental data is available. The only improvement obtained lasing the adaptive grid generation procedure is an increase in the peak of the shock occurring over the cylinder portion of the projectile. However, the peak still does not reach the experimentally predicted value. In an attempt to increase the accuracy in that region, 10 more points were added in the Â£ coordinate direction, thus IM=80, however no further improvement in the solution obtained with the adaptive grid generation technique was evident. The results thus may be the best that the TVD scheme can yield given the current grid configuration and the shortened pressure peak is probably due to the larger dissipative error associated with the TVD scheme. Figure 8.6b shows the grid network corresponding to the converged solution. Again grid point clustering in response to the curvature of the pressure is evident, and as can be seen by the adaptation of the grid network, both shocks occur directly downstream of the expansion waves. TABLE OF CONTENTS ACKNOWLEDGEMENTS ABSTRACT vi CHAPTER I INTRODUCTION 1 II THE ADAPTIVE GRID GENERATION EQUATIONS 7 II.1 The Curvilinear Coordinate System 7 11.2 A Variational Approach 12 11.3 Inversion of the Grid Generation Equations 17 III NUMERICAL SOLUTION OF THE ADAPTIVE GRID GENERATION EQUATIONS 22 III.1 The Newton-Raphson Iterative Method 22 111.2 Modified Solution Procedure 26 111.3 The Adaptive Grid Generation Code 37 TV CONTROL FUNCTIONS 39 IV.1 The Basic Form 39 IV. 2 Smoothing 41 TV. 3 Other Control Functions 42 V FIOW SOLVERS 44 V.l The Governing Equations 44 V.2 Solution Algorithms for the Governing Equations ... 46 V.3 Examples on Fixed Grid Networks 51 VI PRELIMINARY TESTS AND MODIFICATIONS 61 VI. 1 The General Procedure 61 VI.2 Orthogonality 62 VI. 3 Curvature 62 VI.4 Internal Grid Boundaries 69 VII THE SELF-ADAPTIVE COMPUTATIONAL METHOD 76 iv 84 thus, results in too many points being clustered into the boundary layer region at the expense of adequate resolution in other areas. The minimum and maximum grid point spacing along the coordinate lines is 0.00002 and 6.0 respectively. In this direction it is important to have the small spacing at each point along the projectile surface in order to adequately resolve the viscous sublayer region. The parameter a of eq. (4.1.4) was determined for each rÂ¡ coordinate line so that the small spacing would occur for each point on the projectile surface. However, in the Â£ coordinate direction, the shocks and expansion waves are not expected to reach the outer boundary and the small spacing needed along the projectile surface is not required near that boundary. With the minimum and maximum grid point spacing prescribed as 0.02 and 0.22, respectively, a is determined only for the Â£ coordinate line coincident with the projectile surface and that value is used for all the remaining Â£ coordinate lines. This had the effect of increasing the minimum spacing as the magnitude of the shocks and expansion waves diminished near the outer boundary. VIII.2 Results Using the Self-Adaptive Computational Technique: The Beam and Warming Scheme The first application of the self-adaptive computational technique is with the Beam and Warming scheme for the calculation of the projectile flow with sting. The Mach number is 0.96 and the Reynolds number is 76,000, a case for which experimental data is available. The time step used in the calculation is 0.1 and the explicit and implicit dissipation terms are 2 and 4 times the time step respectively. The initial grid network used is that of Figure 5.5 with IM=70 and JM=60. In the first 25 Knowing Axj_j and Ay^j, the new point location can be calculated by . Â£+1 SL+1 solving eqs. (3.1.4) for xj^j and yj^j. This same procedure can be applied to the one dimensional equations for adaptation, eqs. (2.3.9) and (2.3.10). The residue for eq. (2.3.9) is (Rs)i = si+l 2si + si-l + sÂ£pÂ£/p (3.1.9) and the equation for the new point location becomes s{+1=sf+As{ (3.1.10) As{=(Rs)i/2 (3.1.11) The numerical algorithm consists of updating the value of each point on the boundary with eq. (3.1.11) and of sweeping through the domain and updating each interior point using eq. (3.1.7). The criterion used here to indicate convergence is Rmax ^ s (3.1.12) where F^ax is the maximum residue in the interior weighted by the Jacobian Axl,j + Ayf j Rmax= max ( )1/2 (3.1.13) Currently, the criterion 6=0.005 is used. 24 + 1 Z Ax/j xi,rxi,j (3.1.4) 2 + 1 2 AYi/j Yi, jYi,j The derivatives of the residues with respect to the variables and Yi j can be obtained by differentiating eqs. (3.1.1) and are 3(Rl)i,j = -2(ai+a3) , 3xi, j 3(Rj)i,j an,j = -2(b^+b3) (3.1.5) 3(R2),j = -2(C1+C3), 3xi, j 3(R2>,j 3Yi,j -2(d^+d3) An approximate value of AXj^j and AYi,j can be obtained if we require the residue at the next iteration to vanish 2 1 r 3(ri),j 3(Ri),j i/D 3xi, j aYi,j Axi, j 2 3(R2) i, j 3(R2)1,j . Ayi,j . i/j J 3xi/ j aVi, j (3.1.6) Substituting eqs. (3.1.5) into eqs. (3.1.6) and solving for Axj^j and Ay-L, j yields where r Axi,j 1 l *Yi,j J 2 ' (di+d3) -(bi+b3) D . -(C1+C3) (ai+a3) , . (R2)i/j . (3.1.7) D 4((aj+a3)(d^+d3)-(c^+c3)(d^+d3)) (3.1.8) CHAPTER VII THE SELF-ADAPTIVE CDMHJTATTONAL METHOD The adaptive grid generation technique can be incorporated naturally into the solution process of the flow solvers since the flow solver algorithms are time-marching. After a specified number of time steps, the grid network can be adapted to the control functions based on the current values of the flew variables. Thus the time-marching algorithm of the two flow solvers used here, the Thin Layer and TVD codes, proceeds as usual with intermittent calls to the Adaptive Grid code to update the grid network. In this approach, however, the relocation of the grid points in the physical domain due to the grid adaptation must be considered. The values of the flew variables at any time step are defined at the current grid point locations and any movement of the points will therefore distort the solution. This effect is critical for unsteady flow problems since the grid point movement will effect the solution accuracy. For steady flow problems, the grid point movement may effect the convergence rate and if it is large enough, it may cause the solution to diverge. There are two basic approaches to account for the grid point movement resulting when adapting the grid network. The first is to view the grid network as a moving grid and include the time metrics in the coordinate transformation of the Navier-Stokes equations. These metrics can be written in the computational plane as 76 |