UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
NUMERICAL SIMULATION OF THE AERODYNAMICS OF TURBULENT TRANSONIC FLOW PAST A PROJECTILE BY NAEHAUR ERICSHIAU 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 ACKNOWLEDGEMENTS The author wishes to express his sincere gratitude for Dr. ChenChi Hsu, chairman of the committee, for his guidance, support and encouragement during the time of study. Also, his many thanks are extended to the committee members, Dr. Malvern, Dr. Leadon, Dr. Kurzweg and Dr. Shih for their interest and time. The calculation in this study was performed by a Cray system at Bell Laboratories and by a Cray XMP/48 at Pittsburgh Supercomputing Center with all CPU time provided by National Science Foundation through Grant ECS8515761. TABLE OF CONTENTS P ACKNOWLEDGMENTS ........................................ ABSTRACT ............................................... CHAPTER age ii v I INTRODUCTION ................................... II GOVERNING EQUATIONS ............................ NavierStokes Equations ........................ Transformed NavierStokes Equations ............ ThinLayer NavierStokes Equations ............. Axisymmetric ThinLayer NavierStokes Equations .................................... III NUMERICAL ALGORITHM ............................ TimeDifferencing ...................... ....... Approximate Factorization ...................... Numerical Dissipation .......................... Finite Difference Equations .................... Finite Difference Scheme for Axisymmetric ThinLayer NavierStokes Equations ........... IV THINLAYER NAVIERSTOKES CODES ................. Flow Domain .................................... Boundary Conditions ............................. Initial Condition .............................. Viscosity ...... ............................... Forces and Moments ............................. Rate of Convergence ............................ V HYPERBOLIC GRID GENERATION ..................... Governing Equations ............................ Finite Difference Form ......................... VI FLOW PROBLEM ................................... iii VII RESULTS AND DISCUSSION: AXISYMMETRIC FLOWS ..... 59 Inviscid Flow .................................. 59 Flow Domain Determination ..................... 59 Effect of Artificial Dissipations ........... 60 Viscous Flow .................................. 62 Effect of Artificial Dissipations .......... 65 Effect of Eddy Viscosity Distribution ....... 67 Effect of Grid Resolution in Streamwise Direction ......... ........... ............. 67 Effect of Grid Resolution in Normal Direction 71 Importance of Viscous Drag to Pressure Drag 76 VIII RESULTS AND DISCUSSION: THREEDIMENSIONAL FLOWS 89 Flow Cases at 100 Angle of Attack .............. 90 Flow case at M =0.96 ....... ................. 90 Flow case at M =0.91 and 1.20 ...............103 Flow Case at 450 Angle of Attack ...............110 Comparison of Flow Cases at Different Angles of Attack ........................................126 IX CONCLUDING REMARKS ..............................128 APPENDIX ................................................ 130 REFERENCES ................................................ 132 BIOGRAPHICAL SKETCH ..................................... 135 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 NUMERICAL SIMULATION OF THE AERODYNAMICS OF TURBULENT TRANSONIC FLOW PAST A PROJECTILE BY NaeHaur Eric Shiau December 1987 Chairman: ChenChi Hsu Major Department: Engineering Sciences Steady, transonic flow over a projectile was numerically investigated to gain insight and understanding of the thinlayer NavierStokes approximation and to assess the relative importance of viscous effects on the aerodynamic forces acting on a projectile at transonic speeds. The turbulence model used in the computation is an algebraic twolayer eddy viscosity turbulence model proposed by Baldwin and Lomax. The projectile model considered was a secantogivecylinderboattail configuration for which accurate surface pressure measurements are available to assess the numerical results. The sensitivity of the numerical scheme to the grid resolution and added numerical dissipation terms was investigated. The surface shear stress distribution has also been computed to aid in identifying regions of the reversed flow. The computed results, based on a grid system generated by a set of hyperbolic partial differential equations, were verified by comparing with experimental measurements. Good agreement between the numerical results and measured data indicates that the numerical scheme can predict the essential features of the flow field. The calculated results obtained at zero angle of attack and Mach number 0.91 show that the viscous drag is as important as the pressure drag (no base drag). However, the relative importance of viscous drag decreases with increasing Mach number and angle of attack. At 10 and 45 degrees angle of attack, the viscous forces had a negligible effect on both the lift force and the pitching moment. CHAPTER I INTRODUCTION An accurate prediction of the aerodynamic forces is important in the design of flight vehicles. For complex flow fields, an accurate simulation of the flow field is crucial for the successful prediction of such forces. Considerable experimental efforts have been undertaken to enhance our understanding of these complex flows and, recently, with the increase in computational capabilities, there has been a major effort to numerically simulate such flows. The most challenging features of transonic flow to be numerically simulated are the viscous/inviscid and shock/boundary layer interactions. There are two basic approaches used to numerically simulate these flows. One method is a patching technique in which the solution is obtained from inviscid equations and boundary layer equations separately and then patched together in some iterative process. Another approach is to solve a suitable subset of the NavierStokes equations for the entire flow domain. The primary advantage of the patching method is their computational efficiency. However, it cannot solve the separated flow. The latter approach has the advantage that it provides a continuous treatment of the solution through the interaction regions.  1  2  In the transonic flow regime, critical aerodynamic behavior has been observed which can be attributed to the complex shock pattern and local flow expansions that exist over the projectile at these speeds. Shadowgraphs which can be found in reference 7, show the shock pattern existing on a typical projectile configuration for transonic flow. The formation of shock waves imbedded in the flow field produces severe changes in the aerodynamic coefficients such as the drag force and the pitching moment. For example, the nondimensional drag force for a typical projectile shape has been found to change by as much as 100% through a Mach number range of 0.95 to 0.97 [9]. Similar results have also been evident in this study. For such large changes in the magnitude of the aerodynamic forces, it is essential to understand and be capable of computing the features of the flow field that contribute to this phenomenon. In the transonic flow regime, the flow field is characterized by strong viscous/inviscid, shock/boundary layer interactions. It is, therefore, advantageous to use the approach which solves the suitable subset of NavierStokes equations since this method considers these interactions in a fully coupled manner. A three dimensional code was developed in 1978 by Pulliam and Steger at NASA Ames to solve the thinlayer NavierStokes equations for an ideal gas in a transformed boundaryfitted coordinate system [14]. In this code, which is applicable to high speed compressible flow problems, the 3  governing equations are approximated by the spatially factored implicit scheme of Beam and Warming [3,4]. To control the numerical stability of the solution algorithm, second order implicit and fourth order explicit dissipation terms are added. Consequently it results in a block tridiagonal system of equations to be solved at each time level. The turbulence closure model implemented in the code is an algebraic twolayer eddy viscosity model of Baldwin and Lomax [2]. The unsteady thinlayer NavierStokes code has an option for calculating inviscid flow problems and a steady flow solution can be obtained as the converged solution of a time asymptotic flow problem. This code has been widely applied to various flow problems. The threedimensional thinlayer NavierStokes code has also been simplified for axisymmetric flow cases. Both codes have been investigated to some extent by the U.S. Army Ballistic Research Laboratory on transonic projectile aerodynamic problems [9,11]. The planar grid network provided to a NavierStokes code is obtained from a grid generation code GRIDGEN which can give either an elliptic grid or a hyperbolic grid [10]. Moreover, a three dimensional grid system is formed by a sequence of planar grids around the axis of a projectile. For a secantogivecylinderboattail projectile with sting at zero angle of attack, the published result showed that the computed surface pressure coefficient over the secantogive portion and the boattail portion of the projectile agrees 4  rather well with the measured data but the agreement on the cylinder portion of shock/boundary layer interaction region is not very satisfactory for some cases considered. For the projectile model at 2 degrees angle of attack the Cpdistribution computed agrees qualitatively with the measured data but quantitatively the agreement over the cylinder portion and boattail portion is not satisfactory [18]. The published results seem to indicate that the thinlayer NavierStokes codes can give acceptably accurate surface pressure for the complex transonic projectile aerodynamic problem if a good adaptive grid system is provided. However, more precise causes for the unsatisfactory results reported are yet to be investigated; moreover, no result on the skin friction coefficient Cf has been reported and discussed. The GRIDGEN code obtained from the U.S. Army has been investigated and modified [6]. The modified hyperbolic grid is generated using predetermined adaptive boundary points. In order to gain the smoothness property, the hyperbolic grid obtained from GRIDGEN also has been modified to bend the grid lines in the nose region of the projectile, which has been found to be important in the convergence process and for the solution accuracy. The modified hyperbolic grid generation code is efficient, requiring less computer time than elliptic grid generation techniques; moreover, it can provide a good grid network and is therefore used in this study. 5  An accurate prediction of a complex viscous flow field via thinlayer NavierStokes codes depends upon many factors. These include grid resolution, turbulence models, artificial dissipation terms, flow characteristics, etc. For a successful application and/or modification of a NavierStokes code, one must be familiar with certain relevant features of the code. This motivates the present investigation of the computer codes. The main objectives of this study are to further advance our understanding and gain indepth insight on the application of thinlayer NavierStokes codes implemented with the BaldwinLomax turbulence model for accurate computation of transonic projectile aerodynamics. CHAPTER II GOVERNING EQUATIONS NavierStokes Equations In three dimensions, the unsteady compressible NavierStokes equations referenced to Cartesian coordinates without body forces and external heat addition can be written in the following vector form: (1) q 3E 8F aG aEV aFV aG  + + +  + + St ax Sy az ax Sy 8z in which q is an unknown vector, E, F, G are the inviscid flux vectors and EV F GV are the viscous flux vectors. The explicit expressions for these vectors can be written as p pu pv pw E ]T 2 pu pu +p puv puw pv puv pv2+p pvw pw puw pvw pw2+p 0 rxx 'xy 'xz x 0 yx Tyy Tyz y 0 tzx Tzy 'zz z (E+p)u IT (E+p)v ] (E+p)w ]T ]T ]T  6  q = [ E = [ F=[ G=[ EG = [ v= FG= [ G = (2) 7  Here, p is the density; u, v and w are the velocity components in x, y and zdirections, respectively; p is the pressure; and E is the total energy per unit volume which is related to the internal energy e and the kinetic energy by the formula E e + + + (3) E = p [e + (u2 +v2 w2)] (3) For a Newtonian fluid with zero bulk viscosity, the relationships between stresses and velocity gradients are 2 xx 3 2 yy 3 2 TZ  _ zz 3 xy = Tyz = 1 Sx = au 2 ( 2 ax av li (2 ay aw V (2 az au av ay ax av aw az ay aw au ax az aw az au ax av ay (4) )=yx ) = zy ) XZ and the B's in equations (2) are related to viscous dissipation and heat conduction by  8  aT 0x =Uxx + Vxy + xz + k  ax aT y = Uyx + VTyy + Wyz + k (5) ay aT z = UTz + VT + WT + k  asz where p and k are the coefficients of dynamic viscosity and thermal conductivity, respectively, and T is the temperature. The complete set of timedependent NavierStokes equations (1) includes one continuity equation, three momentum equations and one energy equation, but there are seven unknowns p, u, v, w, E, p and T. In order to close the system of equations, the fluid is assumed to be a perfect gas. The equation of state and other relations for a perfect gas are p = pRT e = cT h = T p(6) cv where R is the gas constant, h is enthalpy, and cp and c are the specific heat at constant pressure and volume, respectively. I is the ratio of specific heats, which for air at standard conditions is 1.4. Now, by using equations 9  (3) and (6), one can obtain two additional equations relating p and T to the dependent variables q p = (T 1)[E p (u2 + v + w2)] 2 (7) (71) E 1 T [ (u2 + v2 + w2)] R p 2 Hence, equations (1) and (7) are the governing equations for a class of unsteady compressible flow problems. The number of parameters for flow problems of the same class can be reduced by introducing dimensionless variables. Once the equations are put into nondimensional form, the solution can easily be used to correlate data for succinct presentation using the minimum possible amount of experimental testing. The dimensionless quantities used here are obtained by choosing characteristic scales in the following manner: S x y z t x y z t L L L L/V S u v w T u v w T (8) V V V T E * E p p = V. P.. PV 2 10  where the nondimensional variables are denoted by an asterisk and the freestream conditions by o. For instance, V is a characteristic velocity of the fluid at freestream and L is the characteristic length scale. The substitution of the dimensionless quantities, equation (8) into the governing equations (1) and (7), yields + + + Re ( + + ) (9) at ax 8y 8z ax ay az and 1 *2 *2 *2 p =( 1)[E p (u + v + w)] 2 (10) E 1 2 *2 *2 *2 T = (T1) M [ (u2 + v + w )] p 2 where the Reynolds number Re and the Mach number at freestream are defined as p V L V Re M (11) V'YRT The unknown vector q and flux vectors E F G E , * F G* have the same form of equation (2) as well as the stresses in equations (4) except that all the quantities are nondimensional quantities. The components of P become  11  = * x Uxx * ** * xy "xz * B = U T + V T + W yz zy yx yy yz * S =U t + v T + W z Zx zy zz * +  (Tl)Pr Me 8x * + (T1)Pr M2 y* * + M2 z* 2 * (r1)Pr M az Here, the Prandtl number Pr is cP v Pr p k The typical values of the Prandtl number for air at standard conditions are 0.9 and 0.72 for turbulent and laminar flow, respectively. In order to simplify the notation, the asterisks are removed from the nondimensional equations hereafter except where noted. Transformed NavierStokes Equations To enhance numerical simplicity, accuracy and efficiency, a generalized boundaryfitted curvilinear coordinate system is used. The boundaryfitted coordinate transformation maps all boundary surfaces in the physical domain onto coordinate surfaces of the computational domain. The boundaryfitted curvilinear coordinates not only facilitate implementation of boundary conditions accurately, but also allow all (12) (13) 12  computations to be performed on a fixed square grid in the computational domain. The numerical solution technique, thus, becomes more versatile. Consider a generalized transformation of the form S= ( ( x,y,z,t ) n = n ( x,y,z,t ) 5 = 5 ( x,y,z,t ) = t where (t,n,T,r) are the coordinates of the computational domain and (x,y,z,t) are those of physical domain. It is known, e.g. reference 1, that the Jacobian of the transformation J1, which represent the volume of the grid cells in physical space, is given by jI J xty z + xqyCZt + x;yzn xtyZz xCyRzt (14) x ytzc and the metrics in the computational domain can be related to those in the physical domain by using a chain rule expansion of x,, yt, etc., and then solving for tx, ty' etc., to give  13  tx = J ( y zC yzlI ) ty = J ( zqxC z x1 ) tz = J ( x Iy x Cy ) nx = J ( y cz Ytz ) y = J ( z x zx ) "z = J ( x yt xty ) Cx = J ( ySz ynz ) y = J ( zgx zxt ) Cz = J ( xy, xy ) tt = xSx YjSy z Jz qt = Xqx YXIy zZ z Ct = X Ix YXCy z TC Note that the subscript in the metrics denotes differentiation. The metric identities of the transformation, e.g. reference 1 or 20, are a 1 a St a t a t ()+ ( ) + ( ) +) =0 aa J at J al J at j a t a IX a C x x x ( ) + ( ) +( ) =o at J an J at J (15) a y. a y a y ( ) + ( ) + ( )=0 at J an J at J a 2z a lz a zz z ( )+ ( )+ )= at J an J ac J 14  By applying the generalized transformation to the compressible NavierStokes equations (9) in vector form, dividing by the Jacobian and then rearranging into conservation law form with the use of metric identities, equation (15), the following form of the transformed NavierStokes equations is obtained: aq 8E aF aG aEv F aG + + + = Re ( +  +  ) (16) at as an as as an as in which q = Jq = [ p pu pv pw E ] a8 a8 a8 as E = J ( q + E + F + G ) at ax ay az SJ1 [ pU pUu+Sxp pUv+Syp pUw+Szp (E+p)UStp ]T an an a8 n aa n F = J ( q + E + F + G ) at ax ay az SJ1 [ pV pVu+nxp pVv+n p pVw+nzp (E+p)Vntp T as as ac as G = J ( q + E + F + G ) at ax ay az = j1 [ pW pWu+xp pWv+yp pWW+Cp (E+p)Wetp ]T 0 1 x xx y yx + z zx EV j tx xy + tyyy + tz zy SXX + SyTyz + SzT Zx xz y yz + z zz xOx +tyy z z  15  0 S x xx + yTyx + zTzx Fv Ax xy Ay yyy Az zy x xz y yz A+ zzz x 0 + y y + zBz 0 1 x xx +y yx +z Tzx v xTxy + Cyy z Tzy x Txz + yIyz + z zz Cxx + yy y + zoz The stresses r and temperature gradients in the B terms contain the metrics of the transformation; e.g., au av au 3u au av av av xy=1r()=V(y +y y y x x cx )= xy ay ax at aA ac a& an aS etc. U, V and W are the contravariant velocity components along the three curvilinear coordinates A, 1 and i, respectively. They are defined in terms of the Cartesian velocity components u, v and w as U = tt + U +x + Vy + wz V = At + unx + Vy + wnz (17) W = 4t + ucx + Vy + wz 16  Thin Layer NavierStokes Equations In general, a very large amount of computer time and extensive storage are needed for solving a complex ,aerodynamic problem governed by the complete NavierStokes equations (16). At high Reynolds numbers, however, the flow around the body can be considered to be mostly inviscid, except in a thinlayer region close to the body surface where the viscous effects become important. With this in mind, the thinlayer approximation is used to simplify the NavierStokes equations in generalized coordinates. The thinlayer approximation to the NavierStokes Equations is based on an assumption that the diffusion processes along the body surface are negligibly small in comparison with the diffusion process in the direction normal to the body surface. This is generally true for high Reynolds number flow. As a consequence of this assumption, all viscous terms containing derivatives parallel to the body surface are dropped since they are substantially smaller than viscous terms containing derivatives normal to the surface. Assume that a body surface is mapped into a 4Iplane in the transformed space. Then and T derivatives of the viscous terms can be neglected. Hence, the thinlayer NavierStokes equations in curvilinear coordinates are 17  3q 3E aF 8G 3S 1 + + + = Re (18) at ac aa c where 0 mlu + m2Cx S = 1 mlv + m2 y ml1w + m2Cz mlm3 + m2( xu + yV + z w) S= (x2 + y2 + z2) m2 = P(xuC + yV + 4zw )/3 m3 = (u2 + v2 + w2) /2 + Pr1(l)l(c2 ) and c= J/RT is the local speed of sound. Since the thinlayer NavierStokes equations (18) retain the normal momentum equation, the pressure can vary through the viscous layer. Hence, the thinlayer model is more general than the boundary layer equations and is capable of predicting flow separation. Axisymmetric ThinLayer NavierStokes Equations For an axisymmetric flow in which the flow variables are invariant in the ndirection, the threedimensional thinlayer NavierStokes equations can be further simplified by imposing the following restrictions: (i) the body 18  geometry is axisymmetric; (ii) the flow variables and the contravariant velocities do not vary in the ndirection. Hence, the 8F/ai term of equation (18) can be reduced to a source term H. Based on the relationships between the coordinate systems, the Cartesian coordinates (x,y,z), the cylindrical coordinates (x,O,R) and the curvilinear coordinates (t,n,,) shown in figure 1, one has x = x(.S,C,' ) y = R( t,C,T ) sing z = R( t,(,T ) cos# By evaluating the metric terms with the use of the above relations and substituting them into transformed thinlayer NavierStokes equations (18), the resulting unsteady axisymmetric thinlayer NavierStokes equations become 8q 3E 8G as 1 + + + H = Re (19) at at at 8 With the metrics and flux vectors evaluated at plane 0 = 00, the source vector H [9,11] can then be expressed as  19  Projectiletype body 7Y x = constant plane Li # = constant plane Figure 1. Axisymmetric body and coordinate system  20  H = J1 \ pV[R (UVt) pVR (Vnt) 0 0 + R (WCt)]  p/R 0 which has replaced aF/8n of equation (18). Equation (19) contains only two spatial derivatives but does retain all three momentum equations by allowing nonzero velocity components in the invariant ndirection. Thus, flow over axisymmetric spinning projectiles can be solved with equation (19). CHAPTER III NUMERICAL ALGORITHM There have been various numerical algorithms developed for approximating the NavierStokes equations. For the thinlayer NavierStokes codes obtained from the U.S. Army Ballistic Research Laboratory and used in this study, equations (18) and (19) are approximated by the Beam and Warming scheme [3,4,21], which is an implicit approximate factorization finite difference scheme. This effective finitedifference scheme is described and discussed in the following. Note that for a fixed grid system, the T equals to t. TimeDifferencing Consider three different temporal differencing approximations given by (1) Backward Euler implicit formula n+l n n+1 2 q = At ()l + O[(At)] at  21   22  (2) Trapezoidal formula At aq +q q+l qn [() + ()n+l] + O[(At)3] 2 at at (3) Threepoint backward formula 84 3qn+l 4qn + qn = 2At ()n+ + O[(At)3] at By defining Aqn = qn+ qn, theabove formulas can be combined into a single generalized timedifferencing formula in delta form qn+l 8qn rAt = (1 + s)Aqn sAqn At (1 r) at at + 0[(rs1/2)(At)2 + (At)3] (20) With the appropriate choice of the parameters r and s, equation (20) reproduces many familiar two and threelevel explicit and implicit schemes which are listed in the following table. 23  Table 1. Partial list of schemes contained in equation (20). r s Scheme Truncation error 0 0 Euler, explicit O[At] 0 1/2 Leapfrog, explicit O[(At)2] 1 0 Euler, implicit O[At] 1/2 0 Trapezoidal, implicit O[(At)2] 1 1/2 3pointbackward, implicit 0[(At)2] After substituting equation (18) into equation (20), we have rAt 3E aF Aqn + ( +  l+s at a3i lr aE At ( 1+s 8a aG as aF + + ant as s 1 n+1 n1 Re )n Aq ;c 1+s aG as  Re1 n a4 a which is either first order or second order in time, depending on the choice of r and s as indicated in the above table. The flux vectors E, F and G at the advanced n+1 time step are nonlinear functions of the unknown vector q. It is possible to remove the nonlinear nature and still maintain the order of accuracy of the finite difference approximation by using Taylor series expansion about qn (21)  24  aE En+1 = En + AAqn + O[(At)2] where A =  aq F+1 = Fn + BnAqn + O[(At)2] where B  aq aG Gn+1 = n + CnAq + O[(At)2] where C  aq Matrices A, B and C are so called Jacobian matrices. The viscous flux term can also be linearized by using a method suggested by Steger [16]. In order to apply this linearization method, the coefficients of viscosity p and thermal conductivity k are assumed to be locally independent of q. As a result of this assumption, elements of S have the general form a s = J1 where ai is independent of q = J q given in equation (16), and Bi is a function of q. These elements are linearized in the following manner by Taylor series expansion sin+1 = sin + 1 an [ ( n Aqmn] + O[(At)2 a; m aq 25  These relations can be written in a vector form as Sn+1 = Sn jlnJAqn which has a truncation error of second order in time. The explicit expression of Jacobian matrices A, B, C and M is given in the Appendix. Equation (21) now becomes rAt aA aB aC a [I + ( + + Re(JIMJ)]nAqn l+s at an ac ac s At 3E aF aG as Aq1 ( +  + Re) (22) l+s l+s at an ac a where the matrix I is the identity matrix. Approximate Factorization To increase computational efficiency, approximate factorization is employed for the solution algorithm of equation (22), which reduces the solution process to three onedimensional problems at a given time level. The solution process for Aqn then involves the inversion of three matrices with smaller bandwidth, for which there exist efficient solution algorithms. Consider the following factorized equation approximating equation (22). 26  rAt 8A rAt 3B [ I + [ n +   l+s a l+s an rAt aC a [ I + Re (J_MJ) ] Aqn = RHS(22) (23) l+s ay at Here RHS(22) is used to indicate the righthand side of equation (22). The error introduced by the factorization procedure, equation (23), is the leading third order term At3 aA aB Ba aC aC aA At4 3A aB sC 8qn [( + )+ ]n +O[(At)5] 4 a~ an an ac ac a 8 a3 a3 a8 at which does not affect the second order accuracy of the difference approximation. Numerical Dissipation For a complex aerodynamic problem, a direct application of equation (23) often experiences convergence problems. Hence, dissipation terms are added to the equations to damp out high frequency oscillations and to control numerical instabilities. Several different numerical dissipation terms have been investigated, e.g. reference 13. Currently, the numerical dissipations proposed by Beam and Warming [3,4] and Steger [16] are employed in the thinlayer NavierStokes codes. The second order implicit dissipation operators  27  D = eI j1V A J Dn J1V A (24) D = E J1V A are, respectively, added to the left hand side implicit onedimensional 5, n and C operators of equation (23) while an explicit fourth order dissipation operator DE = E J1[(V A)2 + (VA )2 + (V A)2]J (25) is added to right hand side of equation (23). The symbols A and V are the standard forward and backward difference operators, that is, 1 Afi = (fi+ fi) + O(h) h 1 Vfi = (fi fi) + O(h) h Finite Difference Equations The spatial derivatives appearing in equation (23) have to be represented by finite difference quotients. In the computer codes, central difference approximations are used. Thus, the equation (23) is rewritten as  28  LE LL L Aqn = Rn (26) where rAt L + 6 An + D l+s rAt L I + 6 Bn + DI l+s rAt rAt L I +  6 Cn Re16 (J1nJ) + D 1+s l+s s At Rn Aq (6E + 6 F + 6 G Re1S) + DE l+s l+s At this point, equation (26) can be first or second order accurate in time depending on the choice of parameters r and s. For the finite difference approximation of spatial derivative terms, a second order central differencing 1 6f= (fi fi) + O(h2) 2h is employed on the implicit part of 6 A, 6 B and 6 C for maintaining the block tridiagonal matrix structure, while a fourth order central difference approximation 1 6f (f.+2 + 8fi+1 8fi + f ) + O(h4 12h 29  is used on the convection terms 6 E, 6 F and 6 G of the RHS(26) to keep the convection truncation error from exceeding the magnitude of the truncation error of the viscous terms 6 S which is approximated by second order central difference. In this study, the Euler implicit twolevel scheme is used; i.e., r=l and s=O. Therefore, the solution algorithm is first order accurate in time and second order accurate in space. The unknown vector q at n+l time level governed by equation (26) can now be computed effectively as follows: ** Rn Lq Aq = R ** LE Aq = Aq L Aqn = Aq* qn+l = q + Aqn where Aq and Aq are intermediate solutions. The numerical process thus results in a block tridiagonal matrix inversion for each coordinate. The size of the blocks depends on the number of elements in the unknown vector q. Finite Difference Scheme for Axisymmetric ThinLayer NavierStokes Equations In a manner similar to the threedimensional case, the application of Beam and Warming noniterative implicit scheme 30  to the axisymmetric thinlayer NavierStokes equations (19) results in the following finite difference equations: L Lt Aqn = Qn (27) where s At Qn Aqn1 (6 E + 6 F+ 6 G Re16 S)n 1+s 1+s At Hn E J[(V A)2 + (VA ) n 1+s and Lt and L, are given in equation (26). Again, the Euler implicit scheme, i.e., r=l and s=O, is used in the axisymmetric thinlayer NavierStokes code which results in first order accuracy in time and second order accuracy in space. CHAPTER IV THINLAYER NAVIERSTOKES CODES An axisymmetric thinlayer NavierStokes code and a threedimensional thinlayer NavierStokes code obtained from the U.S. Army Ballistic Research Laboratory (BRL) are employed in this study for the numerical simulation of transonic turbulent flows past a projectile model with a sting. The axisymmetric code is a product of BRL and NASA Ames joint effort [11]; it is a simplified version of the threedimensional code developed at NASA Ames [14]. The application of these codes to transonic projectile aerodynamic problems has been investigated to some extent by the U.S. Army Ballistic Research Laboratory. The threedimensional code and the axisymmetric code are based on the transformed thinlayer NavierStokes equations (18) and (19), respectively, which are valid only for high Reynolds number ideal gas flows. In these codes, the characteristic velocity used in the transformation, equation (8), is the speed of sound at freestream condition c instead of the freestream velocity V ; consequently, the freestream Mach number M in equations (10) and (12) becomes unity. The solution algorithm adopted for the transformed governing,equations (18) and (19) is a noniterative, implicit, approximate factorization, finite difference 31  32  scheme of Beam and Warming which has been discussed in the previous chapter. It is noted that both codes have an option for solving inviscid flow problems while a steady state solution is obtained as the converged solution of the unsteady impulsive flow problem. Both the axisymmetric code and the threedimensional code have been studied in detail. Necessary modifications to the computation of turbulent eddy viscosity have been made to the axisymmetric code and are discussed in the following; furthermore, the code has been vectorized for Cray XMP supercomputers to improve the efficiency of computation. Moreover, essential subroutines for the computation of aerodynamic forces have been written and implemented into both axisymmetric and threedimensional codes. For a successful application and/or modification of a NavierStokes code, one must be familiar with certain relevant features of the code. The essential features of the NavierStokes codes considered and the subroutines implemented are described and discussed in the following subsections for the transonic projectile aerodynamics problems investigated. Flow Domain In this study, flows past an axisymmetric projectile model with a sting are considered for the investigation. For axisymmetric flow problems, a schematic planar flow domain 33  and the corresponding computational domain, ABCD, are shown in figure 2. Here, the line AB is the solid boundary of the projectile model, the line BC is the exit or downstream boundary and the line CD is the outer boundary while the line AD is called the upstream center line. For a projectile model at an angle of attack, that is, a threedimensional flow problem, the flow domain in the physical space is generated by the revolution of the planar domain about the axis of the projectile model. Hence, a grid network for a threedimensional problem can be formed by a sequence of positions of the planar grid ABCD around the axis of the projectile model. The corresponding grid system in the computational domain is then formed by stacking the rectangular grid domain ABCD. Accordingly, as shown in figure 2, the plane ABB'A' is the solid boundary of the projectile model, the plane BCC'B' is the downstream boundary, the plane CDD'C' is the outer boundary, while the plane DAA'D' is a singular plane representing the upstream center line AD. Boundary Conditions For a selected boundaryfitted transformation, such as the one discussed above, suitable boundary conditions must be specified. The boundary conditions in the thinlayer NavierStokes codes are applied explicitly and require the  34  ,., >1 C D *l 0 vC C1 0 00 (0 S.I *4 C *P au >i mI Q 04 35  specification of all variables p, u, v, w and E in unknown vector q on the line ABCD. The application of solid boundary conditions, line AB, is simplified through the boundary fitted curvilinear coordinate transformation since the body surface has been mapped onto the plane of 4=0. The noslip boundary condition for viscous flow is enforced by setting the contravariant velocities to zero, i.e., U=V=W=0 on the projectile surface. For inviscid flow problems, a linear extrapolation of tangency contravariant velocities U and V is used while W=0 for impermeable wall. The velocity components u, v and w are then obtained by inverting equations (17). Densities at the surface are assumed to be equal to the values calculated on the first grid line adjacent to the boundary, which was always well within the viscous sublayer. The pressure along the body surface is obtained from the three transformed inviscid momentum equations. By combining [ xx(&momentum) + Syx(nmomentum) + z x(Cmomentum)] with the use of the continuity equation, metric identities (15) and neglecting viscous terms, one finds (txCx + y y + z z)P + (1x4x + ,yy + +zz)p, + (2 + y2 + z) (28) = P(at + uax + va y + wa ) pU( xU + cyvt + zw) pV(Cxu + yv + CZW) 36  Since viscous terms have negligible effect on surface pressure, the same relation (28) has been used in viscous flow computation [14,16]. Equation (28) is solved at the surface by using second order central differences in the ( and T directions and second order forward differences in the C direction. The total energy E is then obtained from the known pressure, velocities and density at the surface by using equation (10). Along the upstream line of symmetry, line AD, which emanates from the projectile nose, symmetry conditions are imposed by setting the second order finite difference expression for the first derivative equal to zero. In threedimensional problems, the flow variables are determined along the line AD independently from each plane containing line AD. The averaged values obtained from all planes are then used on the singular plane DAA'D'. Along the far field outer boundary, line DC, freestream conditions are imposed. On the downstream boundary, line CB, linear extrapolation is used for all flow variables as the exit plane flow conditions for M 21. For M <1, all flow variables are obtained from linear extrapolation except the total energy E which is obtained from equation (10) by assuming the freestream pressure p. fixed at the exit plane. 37  Initial Condition A steady flow problem solved by the NavierStokes code is treated as an impulsive flow problem. The computational procedure is then started marching in time until an asymptotic steady state is reached. Hence, for the projectile flow problem, the freestream conditions are taken for the initial conditions. To ease the sudden jump in solid boundary conditions when flow impacts the projectile impulsively, the boundary conditions at the surface are turned on slowly over the first 30 time steps. Thus, the boundary conditions of the unknown vector at the body surface are scaled by q = qbc x SL + (1 SL) q4 where SL = (10 15 x r + 6 x r2) r3 r = (NC 1)/30 where NC is the number of iterations, qbc is the correct boundary conditions of the unknown vector which has been discussed in the previous section. 38  Viscosity The effective coefficient of viscosity i in the NavierStokes equations is comprised of two parts, the laminar viscosity py and the turbulent eddy viscosity pt. In the thinlayer NavierStokes codes, the evaluation of laminar viscosity pt as a function of temperature is obtained from Sutherland's theory of viscosity and given by V T T + S 1 +S S = ( )15 = T*5 ( ) O T T + S T +S where the constant S is 198.60R [15] and S =S/T . The turbulent eddy viscosity model employed in the codes is the one proposed by Baldwin and Lomax [2]. It is a twolayer algebraic model in which an eddy viscosity is calculated for an inner and an outer region. The inner region follows the PrandtlVan Driest formulation (Pt)inner = pt21wl (29) with e = ky[l exp(y+/A+)] where 7 is the normal distance to the surface and Iwl is the magnitude of the vorticity w = [(a yaxv)2 + (azvayw)2 + ( zu)21/2 39  and y is the law of the wall coordinate S(Pw w) 1/2 y = Pw where w, P and v are the local shear stress, density and laminar viscosity evaluated at the wall, respectively. The eddy viscosity for the outer region is given by the Clauser formulation (ut)outer = pKCcpFwakeFkleb() (30) where K is the Clauser constant, Ccp is an another constant and Fwake is the smaller value of the following: Fwake = YmaxFmax Fwake = CwkYmaxdif /Fmax The quantities ymax and Fmax are determined from the function F(Y) = ylwl[l exp(y+/A+)] where Fma, is the maximum value of F(y) and yax is the value of y at which it occurs. The Klebanoff intermittency factor Fkleb(Y) is given by 40  Fkleb() = 11 + 5.5(CklebJ7/yax)6 1 and the quantity Udif is the difference between maximum and minimum total speed dif = (u2 + V2 + W2)max1/2 2 + + 2)min/2 where the minimum total speed is zero, unless the spinning projectile is considered. Rather than attempting to match (Vt)inner and (Vt)outer values in the overlap region, both (lt)inner and (Vt)outer are calculated at every grid point in the flow field and the smaller value is used as turbulent eddy viscosity. The constants used for this model are A+ = 26, C = 1.6, Cwk cp Wj 0.25 and Ckleb = 0.3 Note that in this section, all quantities are dimensional except those with a superscript In the original axisymmetric thinlayer NavierStokes code, the turbulent eddy viscosity was calculated only for the first 25 grid points off the projectile surface in order to increase computational efficiency. However, in the current grid networks used, the 25th point is well within the boundary layer and subsequent zero values for the turbulent eddy viscosity in the remaining region of the boundary layer led to instabilities in the numerical algorithm. The code has been modified to include the 41  'turbulent eddy viscosity at every grid point in the entire domain. Forces and Moments The motion of the fluid relative to a solid boundary exerts a dynamic force on the surface which consists of both frictional and pressure forces. These surface forces can create the pitching, yawing and rolling moments relative to the center of gravity. In the cases we considered, the axisymmetric projectile without spinning, there is no contribution to the side (lateral) force, nor the yawing and rolling moments. Note that both pressure and shear stresses can cause the pitching moment for the projectile at an angle of attack. The forces and the pitching moment acting on an arbitrary shape can be obtained by integration of the surface pressure and shear stress distributions and moments about the center of gravity, respectively, as follows: (1) Force components in axial direction FpA = I p sin( dA A (31) FfA = I$ T cos9 dA 42  (2) Force components in normal direction FpN = If p cos8 coso dA (32) FfN = /f T sin8 cos0 dA (3) Pitching moment about the center of gravity MC.G. = II [FN (XC.G. X) + FA Z cos#] dA (33) Here, p and T are the pressure and shear stress distributions on the body. The subscripts p and f refer to the forces due to pressure and shear stresses, respectively, and subscripts A and N indicate the force components in the axial and normal direction, respectively. The angle * starts on the leeward side and measures in the circumferential direction and 0 is the inclination of local slope of the surface, as shown in figure 3. XC.G. is the axial position of the center of gravity measured from the nose of projectile which is 3.6 for this particular projectile. X is the local axial position and Z is local position in the radial direction, measures from the axis of the projectile. The pitching moment is considered positive when it is in the noseup sense. With the same dimensionless quantities in equations (8) and characteristic velocity c forces in the equations (31) and (32) can be put into  43  0 0U >U 'I 0) r C W 0O $4N II U U uU ^! c ,4 ) *C a. , l4 re N 44  Fp = L2 p co F* = F * 1 * Ff = L yi c Ff= Ff Re The ratio of the forces due to pressure and shear stress is * therefore F to Ff /Re. The pressure and shear stress coefficients are expressed as * p p. P p. p ( l / 2 ) p mV 2 ( 1 / 2 ) M .( (35) T /Re f (1/2)p V (1/2)M2 and the pitching moment coefficient as * MC.G. MC.G. Cm (36) (1/2)p V .L (1/2)M. L Because the moment has dimensions of a forcelength product, the additional length scale L in the denominator of moment coefficient is introduced. The component of the resultant force in the direction of the freestream velocity V that passes the body is the drag, which is one of the most important parameters in aerodynamic designs. It can be expressed as  45  FD = FN sina + FA cosa (37) where a is the angle of attack. The corresponding lift force is FL = FN cosa FA sina (38) For axisymmetric flow cases, i.e. a=00, the drag force will be the axial force and lift force will be zero. The computation of forces, moment and shear stress distribution has been added to both thinlayer NavierStokes codes. Thus, the computer codes have become more complete in evaluating the flow features. Rate of Convergence The numerical algorithm is a timemarching scheme in which a steady state solution is obtained as a time asymptotic solution of the unsteady equations. The convergence of the process can be observed from the average residual. (R 2 1/2 Average residual = fm (39) At /(JMAX2)(KMAX2)(LMAX2) 46  where R is either the vector of RHS(26) or RHS(27). The total numbers of grid points in the n and C direction are JMAX, KMAX and LMAX, respectively. And summations over j, k, i range from 2 to JMAX1, KMAX1, LMAX1, respectively. The summation of m is from 1 to the total number elements of unknown vector q, which for this case is 5. For steady state solutions, from equation (26) or (27), we know that the vector R is zero since the solution will no longer vary. This implies then that the average residual becomes zero when the steady state solution is reached. In the numerical process of this study, the convergence criterion for the steady state solutions being reached has been set to the average residual less than 104. CHAPTER V HYPERBOLIC GRID GENERATION In the applications of thinlayer NavierStokes codes, a grid network must be provided. The grid network used in this study is generated from a grid generation code, named GRIDGEN, obtained from the U.S. Army Ballistic Research Laboratory. This code can provide a grid network generated either by elliptic or by hyperbolic equations. For external flow problems, hyperbolic grid generation is feasible. Since specifying the precise positions and shape of the outer boundary is not important, the grid around the geometry of interest can be obtained by the specification of grid positions at one boundary, and then marched outward from this boundary until the far field is reached. The hyperbolic grid generation is a noniterative scheme, so the required computer time is almost equivalent to that of one iterative sweep in solving elliptic grid generation equations. Furthermore, this hyperbolic grid has the orthogonal property, and grid smoothness can be controlled by a selected grid distribution function.  47  48  Governing Equations Construction of hyperbolic planar grid (x,z) proceeds by selecting two differential equations which serve as constraint equations that are solved to obtain the coordinates of grid points [17]. It is known [8] that an orthogonal grid will reduce the truncation error terms associated with finite difference approximation formulated in generalized coordinates; thus one equation selected is that which enforces orthogonality V *0 V4 = 0 or &x For a twodimensional general coordinate transformation, the metric relations and Jacobian J of the transformation are z x z& xt S x z (41) J J J J J = xz4 xcz (42) The equation (40) can then be written in the computational space (&,4) as xx + zz4 = 0 (43) 49  The second equation is chosen to specify a finite grid cell area to insure a nonsingular mapping through the transformation. In numerical implementations, AE=AC=l, the Jacobian in equation (42) is approximately equal to the physical cell area. A set of hyperbolic grid generation equations is then formulated by equations (43) and (42). Equations (43) and (42) are nonlinear equations which can be locally linearized by expanding about a known state (xo,z0). Neglecting the high order term, e.g., (xxo")(xx*),, we have x0x = [x + (x x [ + ( )] = x 4x + xox x x etc. where the superscript 0 denote the known state. Upon substituting these relations into equations (43) and (42), one finds x ox + z oz + x ox + z z = 0 (44) z 0x x oz z0x + x 0z = J + J0 where the grid areas J are to be specified and JO is determined at the known state. Therefore, all quantities of RHS(44) are known. Equations (44) canbe put into the compact form  50  AW + BW = f (45) with x 0 z C B x z f = 0 x X A = B = f= W = z o x 0z o J + 0 z Assume that the C coordinate runs in the direction away from a specified boundary. Then the governing equations for the propagation problem can be written as CW + W = g where C = B1A g = B1f The eigenvalues of C are two distinct real roots which implies that the governing equations (45) form a hyperbolic system. Finite Difference Form In the GRIDGEN code, the finite difference approximations used for the hyperbolic grid equations (45) are second order central difference for t derivatives and first order backward difference approximation for derivatives in the marching direction C; i.e.,  51  ( ),t ( W jlt+ Wjl,t+ ) + o( A2 ) 2A& 1 AC ( WW )1 = ) + O( A ) Then, an implicit finite difference scheme can be formed 1 1 A e W'+1 + + Bj wj,' A2 'W]1l,+1 (46) = ft+ + B, eWi, which forms a 2x2 blocktridiagonal set of equations to be solved on each successive constant Cline. Note that the column matrix f ,e, on RHS(46) is assumed known, since it only involves the quantities J+JO. The grid areas J on the next e+ constant Cline can be specified in any way. In the GRIDGEN code, it has been specified as 1 Jj,+= AL x ASt 2 where AL is the total length of two adjacent cells, as shown in figure 4, and ASe is the specified height of two successive constant 4lines. The distribution AS used in the code is an exponential distribution, and is determined by the relation  52  i:,e I .f Figure 4. Specification of grid area in GRIDGEN code. 53  ASZ = ASO (1 + )i i = 1,IMAX1 (47) where IMAX is the total number of grid points along the given arc length, and ASO is the specified grid spacing at one end of arc. The parameter e is determined by a NewtonRaphson iteration process so that the sum of the above increments matches the known arc length. The grid spacing ASO next to the boundary considered in this study is 0.01 for inviscid flow problems; however, for viscous turbulence cases, first grid spacing ASO of 2x105 is required. Several different distribution functions have been discussed and investigated in detail [19,20,22], and, in comparison, a hyperbolic tangent distribution has the better overall distribution [19,20]. For extensive study of the effectiveness of the grid network on the flow solution, a hyperbolic tangent distribution is implemented in the computer code. A hyperbolic tangent distribution is constructed by 6 1 tanh[( 1) 2 IMAX 1 Si = Stot {1 + tanh( tanh(6/2) (48) ASi = Si. Si = 1, IMAX1 54  where Stot is the total length of the arc and the parameter 6 is determined by 6 AS0 S = (IMAX 1) sinh(6) Stot Two computational planar grids for the projectile configuration with the same boundary distribution and ASO=2x105 grid but different normal distributions, equations (47) and (48), for the projectile configuration are shown in figure 5. We observed that the grid with hyperbolic tangent clustering provides a smoother grid away from the projectile than the one with exponential clustering in the radial Cdirection. However, near the body surface, the grid with exponential clustering is smoother than hyperbolic clustering. A hyperbolic grid generation is a propagation problem; it requires the initial grid points. The grid point distribution on the body surface, in this study, is adapted to the flow solution in the tdirection and then propagated outward in (direction. Furthermore, a three dimensional grid for the axisymmetric projectile geometry considered is generated by rotating a two dimensional planar grid about the axis of symmetry and maintaining a uniform angular space between the rotated planes. 55  2,D I 28 18 8 16e e 18 / (a) Based on exponential clustering. 2/D 18 8 18 18i X' (b) Based on hyperbolic tangent clustering. Figure 5. 90x40 hyperbolic grids with different clusterings in the normal direction. CHAPTER VI FLOW PROBLEM The objectives of this study are to gain further understanding of complex flow simulations using the thinlayer NavierStokes approximation and to investigate the relative importance of the viscous force to the pressure force on a projectile at transonic speed. The geometric configuration considered is an axisymmetric secantogivecylinderboattail (SOCBT) projectile. The projectile model has a 3caliber secantogive part followed by a 2caliber cylinder and 1caliber 7degree boattail as shown in figure 6. The model was supported by a basemounted sting in the experiments. Experimental data are available for the pressure distribution along the projectile surface for a range of Mach numbers in the transonic regime, 0.91, 0.94, 0.96, 0.98, 1.10 and 1.20; at angles of attack of 00, 20, 40, 60 and 100; and at circumferential positions around the model in 22.50 increments [7]. The computational model of the projectile, however, is formed by extending the boattail part of the projectile for another 1.77 calibers to meet the sting. This modification is made to avoid the difficulty of simulating the base flow region. For the external flow problems, which are considered here, the outer boundary needs to be specified. Theoretically, the 56   57  03 C) CO N,, 58  outer boundary should be as far as possible from the projectile, so that freestream conditions can be imposed. The wind tunnel experiments supply the freestream stagnation temperature TO=580R; moreover, the projectile body is assumed adiabatic and impermeable. CHAPTER VII RESULTS AND DISCUSSION: AXISYMMETRIC FLOWS For the purposes of evaluating the accuracy and applicability of the thinlayer NavierStokes codes described previously for complex flow problems, comparisons of measured and predicted aerodynamic characteristics and associated flow fields for various flow conditions have been made. The finite difference scheme for the NavierStokes equations requires the additional artificial second order implicit and fourth order explicit dissipation terms, equations (24) and (25) respectively, to control numerical instabilities. According to a linear stability analysis, the scheme is unconditionally stable for linear equations when EI=2eE [3,4,21]. This weight is used in all subsequent investigation of axisymmetric flow problems. Inviscid Flow Flow Domain Determination For SOCBT projectile model, the experimental data are available in a range of M =0.91 to 1.20. A flow case of M =0.96 is selected for detailed investigation. In order to 59   60  determine adequate flow domains, the axisymmetric NavierStokes equations (19) were solved on a 90x40 hyperbolic grid with exponential clustering, equation (47), in the normal direction. The two domains used, Stot=18 and 24 are about 3 and 4 times the projectilelength, respectively. The results presented in figure 7 show that the computed surface pressures are very similar and give an accurate prediction of the experimental data except on the boattail region. The results indicate that both domains 18 and 24 seem to be sufficiently large for imposing freestream conditions on the outer boundaries. Nevertheless, in theory, the better solution is obtained with a larger domain when enough grid resolution is provided. Effect of Artificial Dissipations Since the flow domain of 18 can provide accurate results, the 78x28 hyperbolic grid with domain 18 was used with two different sets of dissipation terms: one is EI=2E=4At; another is EI=2EE=8At for the inviscid flow case. The influence of artificial dissipation terms on the flow solutions is presented in figure 8. When more dissipation is added, smearing of the numerical solution occurs near the shock/boundary layer interaction' region, about at the middle of cylinder part, since the dissipative nature tends to smear the predicted shock over several mesh intervals. In order to prevent the deterioration of the solution, the  61 M =.96 2. 4. 6. XD Figure 7. Computed surface pressure of inviscid solution with a 90x40 grid but different domain. : Sot=18 ........ ot24 tot'l tot2 M = .96 2. 4. 6. WD Figure 8. Computed surface pressure of inviscid solution with a 78x28 grid but different dissipations.  : EI=2EE=4At : EI=2EE=8At  62  added dissipation terms should be as small as possible; however, the solution algorithm became unstable when CI=2EE=2At is employed. Viscous Flow In turbulent flow applications, the closure turbulence model used in the codes is a twolayer algebraic eddy viscosity turbulence model of Baldwin and Lomax. The computed results obtained show that a 90x60 hyperbolic grid with exponential clustering in Cdirection and flow domain Stot=24, which is about 4 times the projectile length, can give very accurate solutions. Figure 9a shows that the computed surface pressure coefficient Cp is in excellent agreement with measured data while the computed shear stress T, presented in figure 9b indicates that there is no flow separation on the projectile surface. The corresponding pressure and Mach contours are shown in figure 10. The dense regions in the contour plots indicate that large gradients occur; i.e., it implies the corresponding flow property changes rapidly. The coalescence of the Mach lines to the right of the supersonic region represents the position of the shock. Two shock locations, one in the middle of cylinder part and another around the middle of boattail portion, agree quite well with a shadowgraph in reference 7. The expansion waves also are shown at the ogivecylinder and  63  M = .96 2 4 6 X/D (a). Pressure coefficient; .96 0: measured data. 1 2 4 6 X/D (b). Surface shear stress TC Figure 9. Computed results with a 90x60 grid. m w L 0 p (1) ``jV ~27  64  0 2 4 6 (a) Pressure contours p: 0.48+0.86, Ap=0.02 u 2 4 6 (b) Mach contours M: 0.881.0, AM=O.01; M: 1.01.26, AM=0.02 Figure 10. Contour plots of M =0.96. 65  cylinderboattail junctures. The computed nondimensional drag force, D=FDx10 due to surface pressure is Dp=40.28 over the 6caliber projectile, excluding the base drag, while the drag force resulting from skinfriction is Df=20.53. This shows the importance of viscous flow computation for accurate aerodynamic force prediction. Effect of Artificial Dissipations The effect of artificial dissipation terms upon the accuracy of numerical solutions of turbulent flow cases was next considered. The same grid network 90x60 of domain Stot=24 was used but with EI=2EE=8At. The computed C as shown in figure lla seems to agree rather well with the accurate results, while the quality of the solution is degraded near the shock/boundary layer interaction region owing to the nature of dissipation. The resulting pressure drag is Dp=32.91 which is 18% less than that of the accurate solution. Figure llb shows the difference between the calculated shear stress distributions; however, the resulting shear drag is about the same, Df=20.12. The difference in the drag forces due to pressure is mainly due to the difference in the predicted pressure on the boattail part, since the pressure intensity on the cylinder portion does not contribute to the drag force in axisymmetric flow cases. p 66  M = .96 Station 23 39 46 2. 4. 6. oD (a). Pressure coefficient; 0: measured data. M = .96 4 6i iL 2 S2. 4. 6. XoD (b). Surface shear stress C Figure 11. Computed results with a 90x60 grid but different dissipations.  : I=2E=4At : e=2E =8At  67  Effect of Eddy Viscosity Distribution The effect of an averagedtechnique originally implemented in the axisymmetric thinlayer NavierStokes code for computing eddy viscosity also has been investigated on 90x60 grid at M,=0.96. Figure 12 shows the distribution of eddy viscosity with and without the averaging technique at three different boundary point stations identified in figure lla. It is clear that the averaging scheme yields improper zigzag eddy viscosity distribution above the viscous sublayer; nevertheless, it has negligible effect on surface pressure. In fact the computed C distribution is almost exactly the same as that of the accurate solution shown in figure 13; however, the shear stress computed is consistently less than that of the accurate solution over the entire projectile surface. The corresponding pressure drag force is Dp=41.41 and skinfriction drag is Df=18.99, which are different by 3% and 8% respectively. Effect of Grid Resolution in Streamwise Direction The grid networks 78x40 and 90x40 with domain Stot=18 are used for comparison. Of the 12 additional points in the tdirection of the 90x40 grid network, 2 points were added in the ogive part and 10 points were added in cylinder part. From the results presented in figure 14, we observed that the solutions are generally quite the same except in the  68  .. Stat ion M) ...  23 o .. / c ^ "* ~ ' ^ L ..5 ... ..... ......... ...... .  .. ... ... .... S*** .' ... ^'''^^ *^'^'^* 5 8 1 0 1 5 0 Eddy Viscosity (a) Sample eddy viscosity distributions based on an averagedscheme implemented in the original axisymmetric thinlayer NavierStokes code. .15 Station S 23 44 S\ ............ 39 S\  46 .  5o 100 158 ) ,**. / 50 1800 158 Eddy Uiscosity Figure 12. (b) Sample eddy viscosity distributions computed without averaging.  69  M = .96 (a). Pressure coefficient; 0: measured data. I 2. 4. 6. WD (b). Surface shear stress T4 Figure 13. Computed results with and without an averagedscheme. : without averagedscheme : with averagedscheme m a a UI S2 70  Cp M = .96 .4 2. 4. 6. X/D Figure 14. Computed surface pressure with domain S t=18 but different grid resolution in tdirection. ......: 78x40 grid  90x40 grid. 71  cylinder portion which is a shock/boundary layer interaction region. In shock/boundary layer interaction regions, the solutions are more sensitive to grid resolution. The results show that the 90x40 grid, with more points in the cylinder portion, has indeed improved the solution characteristics, especially right after the ogivecylinder juncture. However, the computed Cpdistribution of both cases do not reach the peak value of experimental data on the cylinder part. This could be attributed to insufficient domain or grid resolution on the shock/boundary layer interaction region, etc. Effect of Grid Resolution in Normal Direction To further test for the required grid resolution in the normal direction, the 90x40 grid with domain Stot=24 was solved. The results obtained and shown in figure 15 are compared with the most accurate solution from the 90x60 grid. Both grid networks have the same grid distribution on the body surface and the same exponential clustering function but different numbers of grid points, 40 and 60, in the direction. As is evident in figure 15 by the poor accuracy on the cylinder and boattail portions, the 90x40 grid did not provide enough resolution. The computed drag force due to pressure is 23.95 which is about 41% less than that of the 90x60 case. However, the drag force due to skin friction is 20.28 which has no significant difference from 72  Cp M = .96 .4  8. 2. 4. 6. X/D (a). Pressure coefficient; 0: measured data. M = .96 4 i 2. 4. 6. /D (b). Surface shear stress rT Figure 15. Computed results with Stot=24 but different grid resolution in Cdirection. : 90x40 grid  : 90x60.grid. 73  20.53 of the 90x60 case, even though the shear stress distributions are quite different, as shown in figure 15b. The flow problem with domain Stot=24 has been solved again with a 90x40 hyperbolic grid based on hyperbolic tangent clustering, equations (48), in the normal direction. The characteristics of exponential and hyperbolic tangent clustering functions have been discussed in chapter V; the corresponding grid networks are shown in figure 5. The computed Cp presented in figure 16a shows that the grid resolution using 40 points with exponential clustering in the normal direction is not sufficient; moreover, a grid which resulted from the use of hyperbolic tangent clustering yields a better result; it yielded Dp=27.38 while the other grid resulted in Dp=23.95. These pressure drag forces are, respectively, 32% and 41% less than the accurate value of Dp=40.28. The computed shear stress given in figure 16b indicates that the shear stress distribution is very sensitive to the pressure field, particularly in the shock/boundary layer interaction region; however, the resulting shear drag forces are about the same, 20.05 and 20.28, respectively. It should be pointed out that a grid with hyperbolic tangent clustering does not always give more accurate results. In fact the flow problem with a domain Stot=18 also has been solved on the 78x28 grids. The results obtained and presented in figure 17 show that a grid with hyperbolic tangent clustering gives better result on the boattail portion but the other grid with exponential  74  M = .96 2. 4. 6. XID (a). Pressure coefficient; 0: measured data. 0 *1 2. 4. 6. )WD (b). Surface shear stress T Figure 16. Computed results with a 90x40 grid but different clustering in Cdirection. S : with exponential clustering : with hyperbolic tangent clustering.  75  M = .96 S4. 6. IWD Figure 17. Computed results with a 78x28 grid but different clustering in Cdirection.  : with hyperbolic tangent clustering.  : with exponential clustering  76  clustering yields a much more accurate solution on the cylinder portion of the projectile. Importance of Viscous Drag to Pressure Drag In order to gain further insight as to the relative importance of the viscous drag to the pressure drag, the sequence of transonic flows M =0.91, 0.94, 0.98, 1.10 and 1.20 past the SOCBT projectile model at zero angle of attack have also been solved with the same 90x60 hyperbolic grid using exponential clustering in the normal direction. We observed by comparing computed C distribution with experimental data, shown in figures 18, 19, 20, 21 and 22, that the grid which is a very good grid for M =0.96 apparently is not the one best suited for the other flow cases; however, the results computed are acceptably accurate. The shear stress distributions computed have shown that there is a small region of the reversed flow at the boattail portion for M.=0.91 and 0.94; surprisingly the computed C on the boattail portion is in excellent agreement with measured data. It seems to imply that the BaldwinLomax turbulence model can give an accurate solution for small separated flow regions. The drag forces due to pressure and shear stress are listed in the following table.  77 .91 2 4 6 X/D (a). Pressure coefficient; 0 : measured data. * LJ a 4 822 U) .91 (b). Surface shear stress rTC Figure 18. Computed results with a 90x60 grid.  78 .94 (a). Pressure coefficient; M= 6 C 1 '4 L 0 2 0: measured data. 2 4 6 X/D (b). Surface shear stress 14 Figure 19. Computed results with a 90x60 grid.  79 (a). Pressure coefficient; .98 * : measured data. 6 X/D (b). Surface shear stress rT Figure 20. Computed results with a 90x60 grid. , 4 g2 V)  80 Cp M = 1.1 oo .4 2 4 (a). Pressure coefficient; 0 : measly I M = 1.1 6 W m 4.4 L L o2 u, 0 2 4 (b). Surface shear stress Ir 6 X/D ured data. Figure 21. Computed results with a 90x60 grid. 6 X/D  81 (a). Pressure coefficient; to I) a L U) 0, 0 : measured data. M = 1.2 6 X/D (b). Surface shear stress rT Figure 22. Computed results with a 90x60 grid.  82  Table 2. Drag forces due to pressure and shear stress at different Mach numbers. M 0.91 0.94 0.96 0.98 1.10 1.20 Dp 20.81 32.22 40.28 63.78 78.95 98.24 Df 18.35 19.33 20.53 21.88 27.67 32.56 The data in the table indicate that the forces due to pressure as well as due to skin friction are increased with increasing Mach number. The ratio of the computed pressure drag Dp to skin friction drag Df is, respectively, 1.13, 1.67, 1.96, 2.92, 2.85 and 3.02 for flow cases M =0.91, 0.94, 0.96, 0.98, 1.10 and 1.20. These computed data indicate that the skin friction drag is rather important in the low transonic regime, and is not as significant in the higher transonic regime. The corresponding pressure and Mach contour plots of M =0.91, 0.94, 0.98, 1.10 and 1.20 are shown in figure 23. From a series of comparisons, from figures 10 and 23, we found that initially at M =0.91 only a small disturbance is felt and expansion waves are observed at ogivecylinder and cylinderboattail junctures. Right after the expansion waves, the boundary layer thickness decreases due to the increasing speed and decreasing pressure. As the Mach number increased to 0.94, shocks are seen to occur on the middle of  83  M =0.91 0 2 4 6 (a) Pressure contours p: 0.460.84, Ap=0.02 2 4 6 M: 0.851.0, (b) Mach contours AM=O.O0; M: 1.01.25, Figure 23. Contour plots. AM=0.05  84  M =0.94 2 4 (c) Pressure contours p: 0.460.86, Ap=0.02 Figure 23. (continued) (d) Mach contours M: 0.881.0, AM=O.01; M: 1.01.2, AM=0.05  85  M =0.98 0 2 4 (e) Pressure contours p: 0.460.86, Ap=0.02 8 2 4 6 (f) Mach contours M: 0.881.0, AM=0.01; M: 1.01.26, AM=0.02 Figure 23. (continued)  86  M =1.10 2 4 (g) Pressure contours p: 0.520.92, Ap=0.02 2 4 (h) Mach contours M: 1.01.32, AM=0.02 Figure 23. (continued) 87  M =1.20 2 4 (i) Pressure contours p: 0.5+0.9, Ap=0.02 2 4 (j) Mach contours M: 1.08+1.5, AM=0.02 Figure 23. (continued) 88  the cylinder and boattail parts. Right after the shock waves, the boundary layer thickness increases, which is clearly shown in the Mach contour plot, due to the increasing back pressure and decreasing velocity. With a further increase in Mach number, the regions of supersonic flow have grown and the shocks have moved downstream. At M =0.98 speed, the computed results show that the shock in the middle of the cylinder part has disappeared, and the only shock occurring is at the end of the boattail portion. This is also evident via the shear stress E, plot in figure 20. A further increase in the freestream condition, to M =1.10, shows the beginning of a bow shock at the nose and expansion fans are evident at the surface junctures. At higher speed M =1.20, the contour plots show the shock as well as expansion waves are more inclined toward the flow direction than the case of M =1.10. From the validity of the Baldwin and Lomax turbulence model reported in reference 2, we believe that the above computed results for viscous cases are accurate and good. CHAPTER VIII RESULTS AND DISCUSSION: THREEDIMENSIONAL FLOWS To obtain further understanding on the application of the thinlayer NavierStokes approximation, equation (10), to three dimensional flow problems, the projectile at angles of attack was investigated. The three dimensional grid in this study is obtained by rotating a hyperbolic.planar grid around the axis of the projectile with a uniform angle increment, A0, in the circumferential direction. The angle 0 is the angle measured from the leeward side. It is assumed that the flow problem considered is symmetric with respect to xzplane and thus the grid system generated is only on one half of the projectile. Two extra planes are included to facilitate the calculation of the central difference approximation of the spatial derivatives on the windward and leeward planes. Based on experience in running the axisymmetric code, the 90x60 hyperbolic planar grid from GRIDGEN with domain 24 is used to provide a good grid and ensures obtaining accurate results. This planar grid was rotated to form the three dimensional grid. The total grid points used becomes 90, 19 and 60 in the streamwise circumferential n and radial ( direction, respectively. The transonic flows of M,=0.91,  89  90  0.96, 1.20 past the projectile model at a=100 were chosen for numerical simulation. Flow Cases at 100 Angle of Attack Flow case at M =0.96 For the flow case M,=0.96, the computed results are in excellent agreement with experimental data for the surface pressure distribution. Figures 24 shows the pressure distribution in the longitudinal direction at different circumferential planes. The pressure on the ogive portion is seen to be higher than on the windward side, thus generating a upward force on the projectile. The boattail region, however, shows a reversal with the higher pressure now occurring on the leeward side, thus causing a downward force. The resulting moment about center of gravity causes a noseup behavior of the projectile. This effect can also be seen in the contour plots of figures 25. The computed contour plots also show that the shocks on the leeward side are stronger than those on the windward side. In the Mach contour plot, it is seen that the boundary layer thickness on the leeward side is much thicker than on the windward side, especially after shock which is imbedded in the middle of boattail part of leeward side. The shock positions on the windward side appear to be slightly downstream of those on  91  M= .96 1 0 4  X 7 .  x/ 20 S.. x.*'  X /* / OI ?4 v 2. 6. WO Figure 24. Computed surface pressure at various #stations of SOCBT projectile at a=10.  92  Cp 3 0 "4 ;''\ M = .96 / 3 6 6 x .4. 3   x /6  X d .. 9l   .... .. ."_ .1 A / _0 o 4f .: ^  Figure 24. (continued)  93  M4 00 * II CQ. 0 04 W)a 0* &t * r'  94  0 O + ss 96 .C '._ _. t t 4) 06' 0 rX4 C0 1 O/1 / /// 
Full Text 
xml version 1.0 encoding UTF8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd INGEST IEID EHUCE6PG2_TG0438 INGEST_TIME 20120217T15:56:18Z PACKAGE AA00003797_00001 AGREEMENT_INFO ACCOUNT UF PROJECT UFDC FILES 