Citation |

- Permanent Link:
- http://ufdc.ufl.edu/AA00029914/00001
## Material Information- Title:
- Numerical simulation of transonic turbulent flows over a complex configuration
- Creator:
- Yang, Shu-Cheng, 1954-
- Publication Date:
- 1990
- Language:
- English
- Physical Description:
- vi, 144 leaves : ill. ; 28 cm.
## Subjects- Subjects / Keywords:
- Aircraft wings ( jstor )
Boundary conditions ( jstor ) Grid generation ( jstor ) Jacobians ( jstor ) Navier Stokes equation ( jstor ) Simulations ( jstor ) Solid surfaces ( jstor ) Transonic flow ( jstor ) Turbulent flow ( jstor ) TVD schemes ( jstor ) Aerodynamics, Transonic ( lcsh ) Aerospace Engineering, Mechanics and Engineering Science thesis Ph. D Dissertations, Academic -- Aerospace Engineering, Mechanics and Engineering Science -- UF Turbulence -- Mathematical models ( lcsh ) - Genre:
- bibliography ( marcgt )
non-fiction ( marcgt )
## Notes- Thesis:
- Thesis (Ph. D.)--University of Florida, 1990.
- Bibliography:
- Includes bibliographical references (leaves 139-143).
- General Note:
- Typescript.
- General Note:
- Vita.
- Statement of Responsibility:
- by Shu-Cheng Yang.
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- The University of Florida George A. Smathers Libraries respect the intellectual property rights of others and do not claim any copyright interest in this item. This item may be protected by copyright but is made available here under a claim of fair use (17 U.S.C. Â§107) for non-profit research and educational purposes. Users of this work have responsibility for determining copyright status prior to reusing, publishing or reproducing this item for purposes other than what is allowed by fair use or other copyright exemptions. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder. The Smathers Libraries would like to learn more about this item and invite individuals or organizations to contact the RDS coordinator (ufdissertations@uflib.ufl.edu) with any additional information they can provide.
- Resource Identifier:
- 026233036 ( ALEPH )
25034739 ( OCLC )
## UFDC Membership |

Downloads |

## This item has the following downloads: |

Full Text |

NUMERICAL SIMULATION OF TRANSONIC TURBULENT FLOWS OVER A COMPLEX CONFIGURATION By SHU-CHENG YANG 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 1990 ACKNOWLEDGEMENTS The author wishes to express his sincere appreciation to his dissertation advisor, Dr. Chen-Chi Hsu, for his guidance and encouragement throughout this endeavor. Special thanks are expressed to Dr. W. Shyy for his invaluable criticism of the draft of this dissertation and his stimulating discussions. Many thanks are also extended to the other supervisory committee members, Dr. B.M. Leadon, Dr. U.H. Kurzweg and Dr. A.K. Varma for their interest and comments on the research project. Thanks are also due to Dr. M.H. Chen and Dr. N.H. Shiau for their assistance on their previous work of the TVD scheme and computer codes. Special thanks must go to Mr. J. Ordonez for his tireless proofreading of the manuscript. Also, the author has a debt of gratitude to Dr. N.C. Yao and Dr. C.G. Tu of the Chinese Navy, without whose encouragement and help he would not have had a chance of pursuing the advanced study here at the University of Florida. Finally, the author is indebted to his family for their continuous support. This work was partially supported by a NASA-Ames research grant, NAG2473. Most of the calculations were performed with a Cray-2 of the NAS System at NASA Ames Research Center. ii TABLE OF CONTENTS Page ACKNOWLEDGEMENTS,-.................-.--........-- .-.--....... ii ABSTRACT ... - ---................................................. v CHAPTERS I INTRODUCTION ----......................................... 1 II GOVERNING EQUATIONS ................................... 7 2.1 Navier-Stokes Equations and Equation of State ................ 7 2.2 Non-Dimensionalization -..-................................ 12 2.3 Transformed Navier-Stokes Equations ........................ 14 2.4 Thin-Layer Approximation.--.............................. 17 2.5 Turbulence Model ........................................ 19 III NUMERICAL ALGORITHM - -............................... 25 3.1 Beam and Warming Scheme .--............................. 26 3.2 TVD Scheme .---........................................ 32 3.3 Flow Solver .-- - - --........................................ 43 IV THE FLOW PROBLEM AND SOLUTION STRATEGY .......... 45 4.1 The Flow Problem and Preliminary Work ................... 45 4.2 Solution Strategy .---..................................... 51 4.3 Six-Block Grid Topology - --............................... 54 V GRID GENERATION ----..................................... 62 5.1 Surface Grid Generation ..--............................... 63 5.2 Volume Grid Generation .................................. 66 5.3 The Six-Block Grid ....................................... 74 iii VI SOLUTION METHOD .-.-..-........................-........ 88 6.1 Jacobian and Metric Coefficient ............................. 89 6.2 Boundary Condition .---................................... 95 6.3 Interface Boundary Condition -- --........................... 105 6.4 Convergence Criterion - --.-................................ 108 6.5 Pressure and Shear Stress Coefficient . -..................... 109 VII RESULTS AND DISCUSSION -.-............................. 111 7.1 The Result from Base Grid -............................ 112 7.2 The Result from Refined Grid.-.......................... 115 VIII CONCLUDING REMARKS -- --............................... 130 APPENDICES A JACOBIAN MATRICES .................................... 132 B EIGEN-FUNCTIONS ..--..................................... 134 C CLUSTERING FUNCTIONS ---............................... 136 REFERENCES ................................................ 139 BIOGRAPHICAL SKETCH ...................................... 144 iv 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 TRANSONIC TURBULENT FLOWS OVER A COMPLEX CONFIGURATION By Shu-Cheng Yang December 1990 Chairperson: Chen-Chi Hsu Major Department: Aerospace Engineering, Mechanics and Engineering Science A multi-block computational method is developed for three-dimensional high speed turbulent flows over complex configurations. In this method, the flow domain is divided into several contiguous blocks such that each block is partially bounded by a solid surface. The chosen grid topology has advantages in computing eddy viscosity as well as in applying a thin-layer Navier-Stokes code. In the solution process, each block is then treated as an independent flow problem governed by the same set of time dependent thin-layer Navier-Stokes equations. The interface boundary conditions between blocks are updated at every time step of the solution algorithm. The application of the method with distance weighted interpolation for updating the interface boundary condition is rather simple and straightforward; however, special measures in grid topology and v generation are required. The flow solver employed is a total variation diminishing thin-layer Navier-Stokes code implemented with an algebraic eddy viscosity model of Baldwin-Lomax. In this study, the proposed method is investigated for the simulation of a transonic turbulent flow past a wing-fuselage configuration in which the flow field is properly divided into six blocks. A coarser grid was first used for flow simulation, followed by grid refinement studies. With the use of the refined grid, the calculated wing surface pressure distribution agrees well with measured data. Numerical results obtained also show that the computed shock location and pressure distribution on the wing can be in good agreement with the experiment data if the block grids used are adaptive to important solution characteristics. It is concluded that the proposed method is indeed a very promising approach to be developed further for complex configuration flow simulation. vi CHAPTER I INTRODUCTION For a high speed flight vehicle, numerous complex three-dimensional fluid flow phenomena appear, including the formation and shedding of vortices, the shock-wave boundary-layer interaction and induced flow separation, and the shock on shock interaction. Since the performance of a flight vehicle is strongly affected by these complex flow phenomena, accurate prediction of them and associated aerodynamic forces and moments are critical for advanced design of a modem flight vehicle. At present an accurate wind-tunnel simulation of the complex flow problem in flight conditions can be very difficult, if not impossible. With the advent of supercomputers and the advancement of solution techniques for nonlinear problems, computational fluid dynamics (CFD) is under extensive development and has been applied to complement the wind-tunnel experiment and field test in the design process of a flight vehicle. The use of CFD in the aerodynamic design of flight vehicles has progressed rapidly over the last few years. In the early days, CFD was used to support the validity of a design that developed in the wind-tunnel by trial and error. The state-of-the-art of the CFD method has progressed to the point of being regarded as an important design tool for many flow problems [1]. The reduction of wind-tunnel occupancy is only a direct benefit of CFD. Detailed 1 2 flow information that is unattainable from the wind-tunnel test often can be obtained by numerical simulations. This information about the underlying flows can lead to a better understanding and ultimately to a better design. In fact, the implementation of CFD has fostered a revolution in the design process of flight vehicles. Current CFD methods have only demonstrated an ability to simulate flows about complex geometries with simple physics or about simple geometries with more complex physics. The panel methods are unique in that they provide a capability for solving the flow about completely general configuration. Their principal limitation is that they are restricted to simple physics as modeled by the linear equation. The Euler code provides more flow information than the panel method; however, the important viscous effects are excluded in the approach. For more realistic flow simulations, the application of Navier-Stokes equations is required. However, the use of full Navier-Stokes equations is simply too expensive to be practical for today's computers. An alternative is to use the thinlayer Navier-Stokes equations. In general, the governing equations are approximated by a numerical scheme and then solved in either a structured or an unstructured grid network. The latter is quite flexible in gridding but it requires more work in bookkeeping and computer coding. In the structured grid approach, the governing equations are first approximated by either a finitedifference or a finite-volume scheme, and then solved in a boundary-fitted curvilinear coordinate system. This approach has advantages in boundary condition treatment and computer coding; however, it may have difficulty in 3 generating good grid system for complex configurations. As evidenced by papers published in journals and presented at conferences, e.g. [2-4], the flow simulations involving simple geometries have progressed considerably; however, the capability for turbulent flow problems involving complex 3-D configurations is still in a development state. Hence, further research and development effort is needed to advance the CFD methods for more complex geometry flow simulations. For flow past a complex configuration, such as a flight vehicle, the generation of an acceptable single structured grid for the entire flow field is rather difficult. It is much easier to divide the flow field into several subregions, termed blocks, and then generate a good grid for each subregion. Hence, it is natural to develop a solution method based on the multi-block grid network. There have been a number of multi-block solution techniques proposed for the Euler simulation of aircraft models and other complex configurations [5-7]. However, the topology of block gridding for viscous turbulent flow simulation is more involved than that for inviscid flows and requires a lot more grid points to resolve the thin viscous layer. In general, a good block grid topology for Euler simulation may not be a proper one for Navier-Stokes simulation. There are also zonal methods developed to use different equation sets in different regions of the flow field. In general, the flow field is divided into inviscid and viscous zones, and then the Euler equations are applied in the inviscid zone, while the boundary-layer or Navier-Stokes equations are used in the viscous zone. A zonal approach that was based on the Euler/thin-layer NavierStokes equations was proposed and tested successfully for transonic wing flow [8]. 4 Another zonal algorithm that applied the Euler and thin-layer Navier-Stokes equations for simulating the flow field of isolated wings was reported in [9]. The flow field was divided into four zones; two inviscid zones with coarse grids and two viscous zones with clustered grids. The computed results were in good agreement with experimental data for the surface pressure, except in the immediate vicinity of the tip and in the shock-induced separated region. The zonal approach also has been used to simulate flow problems for more complex geometry. A transonic viscous flow past an F-16A wing-fuselage configuration has been simulated by a zonal approach [10]. The flow field was divided into as many as 16 zones; in the inner zones adjacent to no-slip surfaces the thin-layer Navier-Stokes equations are solved, while in the outer zones the Euler equations are used. The prediction of wing surface pressures was quite good, except at the leading edge and the aft-shock position. The zonal approach has been shown to be rather efficient. However, the difficulties are to locate properly the interfaces for the inviscid and viscous zones and to patch the solutions that are obtained from different equation sets at these interfaces. Recently, a zonal method with Chimera overset grid scheme was investigated for subsonic turbulent flow about the F-18 fuselage forebody and the combined wing-fuselage [11]. Special coding was required for the overlapped grid region, but the computed results have been shown to be in good agreement with flight-test flow visualization and surface pressure measurements. In this study, a novel multi-block computational method is proposed and explored for the simulation of transonic turbulent flows over complex 5 configurations. In this method, the flow domain is divided into several contiguous blocks such that each block is partially bounded by a solid surface. For ease in applying thin-layer approximation and the Baldwin-Lomax turbulence model, the solid surface is mapped onto an entire boundary plane in computational space. In the solution process, each block is then treated as an independent flow problem governed by the same set of time dependent thin-layer Navier-Stokes equations. The interface boundary conditions between blocks are updated at every time step of the solution algorithm. The solution algorithm with distance weighted interpolation for updating the interface boundary condition and the computer coding are rather simple and straightforward. The multi-block method offers several distinct advantages over the single block computation. In fact, the use of block grids makes the problem of simulating flow fields about complex geometry more tenable. If a geometric feature is changed, only the related block requires to be modified without completely redoing the basic grid topology. Also, the solution procedures are adapted to the modern computer architecture; in parallel processing, each block can be solved at different processors; or in case short of main memory, each block can be solved at a time while the other blocks remain on extended storage. A transonic turbulent flow over a wing-fuselage configuration is investigated for the development of the proposed method. The flow field is properly divided into six contiguous blocks. A coarse base grid was first generated and used for the flow simulation, and then a refined grid was generated from the base grid by halving the grid spacing along the wing in the 6 chordwise direction. A flow simulation package that includes grid generation, flow solver, and plotting program, was developed and applied to the six-block flow problem. The method was first investigated for a transonic turbulent flow of Mach 0.8 past the wing-fuselage configuration at zero angle of attack, then simulations for different angles of attack, a=4* and a=-3', were conducted to gain further insight and understanding on the solution algorithm. With the use of the refined grid, the calculated wing surface pressure distribution agrees well with measured data. Numerical experiments conducted also show that special measures and care are required in generating good grid systems involving excessive distortion between the physical and computational domains; however, the results obtained indicated that the proposed method is a very promising approach for the simulation of complex configuration aerodynamics. CHAPTER II GOVERNING EQUATIONS Some theoretical background that required for transonic turbulent flow simulations are described in this chapter. The Navier-Stokes equations and ideal gas equation of state are first discussed. Then, they are non-dimensionlized and transformed to a curvilinear coordinate system for ease in numerical applications. The thin-layer Navier-Stokes equations valid for high Reynolds number flows are obtained by neglecting the diffusion terms along the direction of the solid surface. The time-averaged thin-layer Navier-stokes equations are used for turbulent flow simulations and the Baldwin-Lomax model [12] is used as the turbulence closure. 2.1 Navier-Stokes Equations and Equation of State 2.1.1 Navier-Stokes Equations The fundamental equations of fluid dynamics are based on the following universal laws of conservation: (a) conservation of mass; (b) conservation of momentum; (c) conservation of energy. The Navier-Stokes equations that govern most of the compressible flow problems can be obtained by applying those universal laws to a fluid flow [13-14]. 7 8 Continuity equation. The conservation of mass law applied to a fluid passing through an infinitesimal, fixed control volume without sources or sinks of mass yields 8p + apu + 8pV + =pW 0 (2.1) at ax ay az where p is the fluid density and u, v, w are the components of fluid velocity in x, y, z direction, respectively. Momentum equations. The conservation of momentum law applied to the control volume while neglecting body forces reads apu+ (PU2+p)+ apuv apuw _ a. ar a r at ax ay az ax ay az apv + apvu + a( apvw ar ar ar - y - + -" y+ - I-z (2.2) at ax ay az ax ay az apw apwu + apwv + a _ ar, +r ar. at ax ay az ax ay az where the components of the viscous stress tensor r, for a Newtonian fluid are given by the constitutive equation au a w= 2 au av aw Tx1C = ax ay a ax 3 A(-+ -- + - )+2 - (2 av a) ax az ax 3 ax ay az av aw au av 2 a w a r=A(--+--+--)+2,u - =-(- ---ay az ax ay 3 ay az ax aw au av aw = 2 aw au av r= X( -+ -+--)+2 - -p( - --) (2.2a) az ax ay az 3 az ax ay au av r= 1( -+ - ) = r S ay ax 9 av aw 8z ay aw au rz = p( + - ) = *r, ax az Here, p is the thermodynamic pressure and y is the dynamic viscosity. The fluid is assumed to be isotropic and satisfying the Stokes condition, A = -2/3 p, which is a good approximation for flow problems without relaxation process. For air, the viscosity p depends mainly upon the temperature and can be estimated by using an interpolation formula based on Sutherland's theory of viscosity, i.e., T m T. + S = T + S (2.2b) where p. denotes the viscosity at the reference temperature T., and S is a constant which assumes the value 198.6'R [15]. Energy equation. The conservation of energy law applied to the control volume and assuming no external heat addition results aE a a a = #3 ~ 23 + [(E+p)u]+ [(E +p)v] + [(E+p)w] - Ip" +g -L- (2.3) at ax ay az ax ay az where aT , = Urx +vr +wr +k -ax aT + = uryx+vryy+wry +k ay (2.3a) ay aT 16Z = urz +vrT,+wr. +k -ay 10 Here, E is the total energy per unit volume, T is the temperature and k is the thermal conductivity. The heat loss by conduction qj is related to the temperature gradient by Fourier's law of conduction aT -k -T (2.4) axj The equations described above, Eqs.(2.1-2.3), constitute a complete set of the time dependent Navier-Stokes equations that is valid for a class of compressible flow problems without heat addition and chemical reaction (or relaxation). The Navier-Stokes equations includes one continuity equation, three momentum equations and one energy equation, but there are seven unknowns, namely, p, u, v, w, p, E, and T. The relationships between the thermodynamic variables (p, p, T, E) are established by the equation of state of the working fluid. 2.1.2 Equation of State In this study, air is the assumed working fluid and the perfect gas equation of state is employed. For a perfect gas, the thermal and caloric equations of state are p = pRT, e = CT ( C, = constant) (2.5) respectively. Here, R is the universal gas constant, e is the internal energy. Some other useful relations are h = CT, -y = C,/C,, a2 = yRT (2.6) where h is the enthalpy, C, and C, are the specific heat at constant pressure and volume, respectively, y is the ratio of specific heats, and a is the local speed of sound. 11 If only internal energy e and kinetic energy are considered significant, the total energy E can be written as E = p[e+(u2+v2+w2)/2] Consequently, one obtains two specific equations relating p, T to the dependent variables p, u, v, w, and E, i.e., p = (-y-1)[E-g 1P(U2+V2+w2)] p 2 [ (2.7) T = _[ Ef..12+V2+wlv1 R P 2 If expresses heat conduction q in terms of internal energy through the use of caloric equation of state, one can rewrite Eq.(2.4) as follows. aT k ae yp ae q=-k -- - (2.8) ax C, axj Pr axj Here, the Prandtl number Pr is defined as Pr = Cp/k, for air Pr = 0.72, approximately, and the internal energy is given by E 12 e = P- y(u +v2+w2) (2.8a) and the p's in Eq.(2.3a) can also be rewritten as fx = urx +vr, +wr, + y ae Pr ax fi = ur +vr+wr, + / (2.9) '~Pr ay fi = ur +vr3,+wr, + -Y ae Pr az 12 2.1.3 Vector Form The Navier-Stokes equations, Eqs.(2.1-2.3), form a system of five equations and can be cast in vector form as aq aE aF 8G aE, aF, aG, -_+ -+ -+ -- - + --+ - (2.10) at ax ay az ax ay az where q is an unknown vector, E, F, G are the flux vectors and E,, F,, G, are the viscous flux vectors. The explicit expressions for these vectors are pu P p2 +p pvU q pv, E= puv F pv2+p ,E (E+p)u J (E+p)v (2. 10a) pw ' O ~~ '' pwu r G = pwv - E,= rX, F, = rYY- G, = y G1pwv+ rx Ez z (E + p~w. lax 1,. A in which ru are given by Eq.(2.2a), p's are defined in Eq.(2.3a), and p, T are related to q by Eq.(2.7). 2.2 Non-dimensionalization To reduce the number of parameters that characterize the same class of problems, a dimensional analysis is employed. The non-dimensional quantities are obtained by choosing characteristic scales in the following manner: . x . y * = * t L ' ' L ' L/a, . U . v . w T U= -, V - , w T = T - (2.11) a. a. a.T. 13 E . p . p E= P= -, p = 2 where the non-dimensional variables are denoted by an asterisk *, the freestream conditions by -, a, is the freestream speed of sound, and L is a characteristic length scale. By substituting the above relations into the governing equations, Eq.(2.10), the non-dimensional equations in vector form are aq aE aF* aG* aE aF'* aG -- + --. + - + -r = Re-'( -+ -+ -i) (2.12) at ax ay az ax ay az where the Reynolds number, Re, is defined as Re = and the relationships from the perfect gas equation of state, Eq.(2.7) become 1 * p (--1)[E*-,P(u +v"+w)= (_y-1)p e' E* 2 T = 7(7-1)[.f* 2.(u2+v2+w")] (2.13) The unknown vector q* and flux vectors E*, F*, G*, E,*, F,' and G, have the same form as in Eq.(2.10a), except that all the quantities are now non-dimensional. In summary, the dimensionless equations govern a class of flow problem. Only two parameters, namely, the Reynolds number, Re, and the Prandtl number, Pr, appear in the equations, the Mach number is dropped out due to the use of a, as characteristic velocity scale. In the following, the asterisks * are removed from the non-dimensional equations except where noted. 14 2.3 Transformed Navier-Stokes Equations In the finite-difference approach, to enhance numerical accuracy and efficiency in applying boundary conditions, the governing equations are transformed from the Cartesian coordinates to a boundary-fitted curvilinear coordinate system. Consequently, a general computational code can be developed such that the computation is performed in an uniformly-discretized computational domain [16]. If one introduces a generalized coordinate system = ( xy,z,t ) 17 = ( x,y,z,t ) (2.14) = ( x,y,z,t ) r =t and applies the chain rule of partial differentiation, the non-dimensional governing equation, Eq.(2.12), becomes q,+ q,+ q,,, + q,+ (,Ec + ,E,7+ .Ec) + ( F+ +yF,, ++F) + (Gc+,G,,+ G ) = Re"[( .E,+,E,,+ E,,)+( yF,,+ +yF,,+ F,,)+( G, ++,G,,+,G )] (2.15) The metric coefficients appearing in these equations can be determined by the matrix relation as follow [C '7Y ?z7t yf y,, yC y' (2.16) x~~ z z z, 0' y7z7 = 1 dy y 1y1y' or 15 [( 1 x x,, x C 17, 1, = J, y y . f y . I Iz f z,, z C [y,,ze-yrz,, -(x,,ze-x.z,,) x,,y,-x~y, -( z-yze) xzz, -x-z -(x"yC-x-yX) yez,,yze -(x z,,-x,,7z ) x~y,,-x,,y ctl (xt ( 17 ~ f 1 y7 , 17 7 where J is Jacobian of the transformation and J1, representing the volume of the grid cell in physical space, is given by 1= a(x,y,z) x x, , x a(Q,1S) I zc z, z' = x +yz-y z,)-x,(yz -yrzd)+ x(yez,-y,,ze) (2.17) In order to put the transformed equations, Eq.(2.15), into conservation-law form, the equation is first divided by the Jacobian and then rearranged to give (q,)/J + (q + Ee + Fc+ zGe)/J + ( +,+q,, ++,,+ 7F,+ zG,)/J + (q&+ .Ec+ YFc+ zG )/J = Re-[((E,,+ X,+ezG,,) + ( + +E,+tF,+%G,) + ( .E,+ YF,+ zG,)]/J and then rewritten into the following conservation law form (q/J), + ((q+ .E+ YF+ zG)/J)c + ((,7q+t E+yF+nzG)/J), + (( q+ .E+ YF+ zG)/J)c = Re-'[((QE,+ YF,+ zG,)/J)c + + ((7E,+ ,7F,+ 9,G,)/J),, + ((,E,+ F,+ zG,)/J)j] (2.18) in which the following metric identities [10] have been applied. ()+ ( )+ ( ) = 0 ( )+ ( )+0 a ~~ ~ 0 , J a 16 a J (-1 )+ an J (2.19) a ( ) = 0 a J (L)+ a(-)+ a ( I.) = 0 a8 J an J a J The transformed equation, Eq.(2.18), has the same form as the original equation, Eq.(2.12), and can be written as AA A A A aq +E aF aG . E, -+-+ - + - = Re(--+ 8r a8 an a a A aF, + - + an aG, --) (2.20) where A q = q/J, F = (ji~q+q.E+ryF+rG)/J, A E, (c.E,+ CYF,+ .G,)/J, A E = (,q+ ,E+ YF+ ,G)/J, A G = ( tq+ ,E+ ,F+ ,G)/J, A F, = (q.E,+P,~.F,+n.iG,)/J, G, = ( .E,+ gF,+ G,)/J The explicit expressions for the transformed unknown vector, flux vectors and viscous flux vectors are A pu q =J-' pv , pw I E PV ) PuV + 17,P F = J-P VV+17p -, (E+p)V-ri7p . 0 X X+ y X+ rzx C.Pr + IYY z)% pU E = J pvuU+C p pwU+ p (E+p)U-QpJ AP pUW + = *1puW+gp G pvW+ryp -y, pwW+ p (E+p)W- tp I , A F, [00 =J-1 IlxTXX + 7yryx + 7zrzx 1 7x7rxy +17yryy + nz1zy ??xr z+ 7 yT + nzrT ,7'x +VV +770 (2.20a) (2.20b) I , A V F., 17 0 AX i %D.Y Yryy+ *zzi =Xp J1+ UY ,,+ ZA where U, V and W are contravariant velocity components given by U = 1+c'u+q+ 'w V = ,+17,u+iyv+,7w (2.20c) W = r+ Xu+ Yv+ w It is noted that the velocity gradients in r,, Eq.(2.2a), and the temperature gradient in p's, Eq.(2.3a), have to be transformed into computational space also, e.g., 'rx au +av +8w )+A u r= -+ + -)+2pax ay az ax au u au +av _av _v a8 a a a Y an a + ( aw-,L + aw nz+ aw .)]+2u[( a L + a ,+ -au at an a a an as 2.4 Thin-Layer Approximation The numerical simulation of complete Navier-Stokes equations, Eq.(2.20), is in general rather time consuming and needs a huge amount of computer storage. For high Reynolds number flows, however, the viscous effects are often confined to a thin layer near no-slip boundaries, called the boundary-layer, outside of which the flow can be considered to be mostly inviscid. With this in mind, it is reasonable to assume that the diffusion processes along the body 18 surface are negligibly small in comparison with those in the direction normal to the body surface. 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, while all other terms in the momentum equations are retained. One of the advantages of retaining the terms which are normally neglected in boundarylayer theory is that separated and reverse flow regions can also be computed. Also, flows which contain a large normal pressure gradient can be readily computed. Another justification of the use of a thin-layer approximation can be given from a computational viewpoint. In practice, a substantial fraction of the available computer storage and time is expended in resolving the cross-stream gradients in the boundary layer since a highly stretched grid is required, and the resolution along the body is treated as what is needed in inviscid flow. As a result, even though the full derivatives are retained in the equations, the gradients along the body are not resolved in an adequate manner. Hence, one can drop those terms that are not being adequately resolved, provided that they are reasonably small. Since the thin-layer model requires a boundary-layer type of coordinate system, it is assumed that the solid body surface is mapped onto an entire g1plane in the transformed space, then the thin-layer approximation requires that the and t derivatives in the viscous terms be dropped. Upon simplifying the complete Navier-Stokes equations, Eq.(2.20), the system of governing equations becomes 19 aq a E aF aG .1 S + - + -+- = Re( --) (2.21) ar a a7 a a where 0 J ~ m~u,+ m27. S =J miv + m2y- (2.21a) m1w. +W.2~ I = A( .2+ ,2+ .2) m2 = u(,u, + Yv,+ zw )/3 (2.21b) m3 (u2+v2+w2),/2 + Pr-1(-11)-(T), Note that the notation ^ has been removed from the above equations for simplicity and the heat conduction term are replaced by q _ y ae _ aT (2.22) Pr a (-y-1)Pr a8 The thin-layer Navier-Stokes equations, Eq.(2.21), in conjunction with the equation of state, Eq.(2.13), constitute the required governing equations for this study. They are valid for both laminar and turbulent flow problems. However, due to the widely disparate scales of turbulence, it is impractical to apply the equations directly for turbulent flow simulation. An alternative is to use the time-averaged equations coupled with a turbulence model. 2.5 Turbulence Model The time-averaged equation is obtained by decomposing the dependent variables into time mean and fluctuating components and then time-averaging the 20 entire equation. In the averaging process, new unknowns will be introduced and must be treated by some approximating techniques so as to close the problem. In this study, the mass-weighted time-averaged equations are used for turbulence flow simulation with the eddy viscosity that appears in the averaged equations provided by the Baldwin-Lomax two-layer algebraic model. In this approch, the mass-weighted time-averaged thin-layer Navier-Stokes equations will have the same form as the original equations, Eq.(2.21), while all the dependent variables are time-averaged. The effective coefficient of viscosity I is comprised of two parts, the laminar viscosity p, and the turbulent eddy viscosity p,. The laminar viscosity p, is obtained from Sutherland's theory of viscosity, Eq.(2.2b), and the eddy viscosity p, is provided by a turbulence model [16]. Similarly, the effective thermal conductivity k is also comprised of two parts, k, and k,. The heat conduction is expressed in dimensional quantities as aT C, & a q = -(k,+k,) - _( C,+ ) (2.23) az Pr, Pr, az The heat conduction term in Eq.(2.22) is then replaced by its dimensionless counterpart q - ( A' + At ) aT (2.23a) q - 1 Pr, Pr, ag Here, the turbulent Prandtl number Pr, is defined as Pr, = pcp/k, and is about 0.9 for air [10]. In this approach, the eddy viscosity pt is the only quantity that has to be modeled. In what follows, the theoretical background about the algebraic model will be discussed first, and then the Baldwin Lomax model is described. 21 2.5.1 Two layer algebraic model In analogy with the coefficient of viscosity in Stoke's law for laminar flow, Boussinesq [17] introduced a mixing coefficient, pt, or referred to as eddyviscosity, for the Reynolds-stress that appears in the time-averaged equation. au t= a (2.24) ay where u is the streamwise velocity, y is the spatial coordinate normal to the solid surface, and au/ay is assumed to be the only significant component of strain rate. With Prandtl's mixing length hypothesis, the coefficient pi, is then related to mean velocity field by p = P12 -u (2.25) ay where 1 is the mixing length. The algebraic model is then to establish a relation between the mixing length and the characteristic length of the respective flow. For simple wall shear flows, the mixing length in the law of the wall region is proportional to the normal distance from the wall, i.e., 1 = ky, where k = 0.41 is a universal constant determined experimentally. In the viscous sublayer, however, the mixing length theory is not valid. A corrected form of the mixing length valid down to the wall was deduced by van Driest. He introduced a damping factor and wrote 1 = ky[1-exp(-y+/A+)] (2.26) 22 where y' = ypu,/ is the law of the wall coordinate, A' denotes an empirical constant which equals 26 for flat-plate flow, and u, the friction velocity, is defined as u, = (r./p)', where r, is the shear stress on the wail. In the law of the wake region, the Clauser model is commonly used in conjunction with a factor to account for the intermittent effect. The eddy viscosity is then written as y, = kCUs*Y (2.27) where the value k, is obtained experimentally and assumed to be 0.0168, U is the streamwise velocity on the edge of the boundary layer, 6* is the boundary layer thickness, -y is the intermittency factor deduced from a curve fit of measured data and given by Reynolds and Cebeci [18] as y = [1+5.5(y/6)6]-l (2.27a) where 6 is the boundary layer thickness. Based on the above arguments, Cebeci and Smith [19] proposed a well known two-layer algebraic model as follows. In the inner-layer: t= p12 au (2.28) ay 1 = ky[1-exp(-y+/A*)] (2.28a) In the outer-layer: p= kCU6*-y (2.28b) -y = [1+5.5(y/6)6]1 (2.28c) 23 This model has been shown to be in excellent agreement with experimental observations for simple shear flows. However, the model requires the boundary layer thickness and the edge velocity which constitute a practical disadvantage in the implementation of the model. In an attempt to overcome this difficulty, Baldwin and Lomax [12] modified the outer layer formulation to eliminate the necessity of finding the boundary layer thickness. In addition, they replaced the velocity gradient by the vorticity, which is invariant under coordinate transformation, so that the model can be easily applied to 3-D problems. The Baldwin-Lomax model is described next. 2.5.2 Baldwin-Lomax Turbulence model The Baldwin-Lomax model is similar to the Cebeci-Smith two-layer algebraic turbulence model. The eddy viscosity in a profile is calculated for each inner and outer layer, and then matching both layer as follows [12]: t I(Ji..,, , y5 , (2.29) ("I)outer , y z y, where y. is the smallest value of y at which values the inner and outer formulas are equal. The inner layer follows the Prandtl-Van Driest formulation ('t)ir.', = pl21wI with 1 = ky[1-exp(-y*/A')] (2.30) where y is the normal distance to the surface and Iw is the magnitude of the vorticity Jw| = [(u,-v.)2+(v.-w,)2+(w.-u)2 ]' (2.31) 24 y+ = y(pr.)'/p, is the law of the wall coordinate, A+ = 26, and p., r., A, are the local density, shear stress, and laminar viscosity evaluated at wall, respectively. The eddy viscosity for the outer layer is given by (pt),.t,, = pk CF.IkFkib(Y) (2.32) where k. = 0.0168 is the Clauser constant and C, = 1.6. The parameter F., is given by F~k = min { F"m / where Ck = 0.25 (2.32a) CAgy.U'djf/F. and Udi = (u2+v2+W2)1,. - (U2+V2+W2)",,i. (2.32b) The quantities y. and F, are determined from the function F(y) = ylwl[1-exp(-y+/A')] (2.32c) where F. is the maximum value of F(y) and y. is the value of y at which F. occurs. The Klebanoff intermittency factor Felb(y) is given by Fkleb(y) = [1+5.5(Cjy/yM)']4 and Ckeb = 0.3 (2.32d) The Baldwin-Lomax model has been investigated and applied to a variety of 2-D and 3-D flow calculations with satisfactory results [10]. The model is capable of handling the compressible flow problems with modest flow separation [20]. The implementation of algebraic models requires only a small amount of computer time and storage; this is particularly important in 3-D computations. CHAPTER III NUMERICAL ALGORITHM In numerical simulations, the transformed governing equations described in the previous chapter are approximated by finite-difference equations and then solved in a equally spaced computational space. The basic computer code used in this study is a TVD thin-layer Navier-Stokes code which is based on an original version of thin-layer Navier-Stokes code, ARC3D of NASA Ames Research Center. The original ARC3D was based on the Beam and Warming scheme [21], and has been improved for its efficiency, accuracy and convergence [22]; for instance, there are options for using different numerical dissipations for better stability criteria. Alternatively, to improve the robustness of the code, a symmetric Total Variation Diminishing (TVD) scheme was implemented into the original ARC3D. This modified ARC3D has shown to be more efficient and robust than the original one [23-26], and hence is adopted in this study. The Beam-Warming scheme and the TVD scheme were constructed by totally different philosophies; however, in applications they can be viewed as the same scheme with different numerical dissipation terms [27]. Both schemes will be described and discussed in the following. 25 26 3.1 Beam and Warming Scheme The Beam-Warming scheme is an non-iterative, implicit approximate factorization finite-difference scheme. The implicit operator is factorized (spatially split) to obtain an ADI-like ( Alternating Direction Implicit ) algorithm such that the three-dimensional inversion process is reduced to three onedimensional processes. The basic algorithm is second-order accurate in both time and space [28]. The implementation of this scheme to the system of governing equations, Eq.(2.21), is described as follows. 3.1.1 Temporal-Differencing To approximate the time derivative of Eq.(2.21), the well-known secondorder Trapezoidal Rule or Crank-Nicolson Formula is employed to give Aq" 1 aq"*l 8q" - = ( + q) + O[(Ar)2] where Aq'=qn+l-qn (3.1) Ar 2 ar a r If one rewrites Eq.(2.21) as aq aE aF aG - aS ar a ai a as and then substitutes it into Eq.(3.1), the following equation results Aq r aE aF .kIGRe(aS ,q" - { --- + - + -Re-'( --)]"+1 2 a an a a aE aF aG 1aS + [ + + --Re-'( )]"} + O[(Ar)3] (3.2) a a?? a a In the above equation, the flux vectors E, F, G, and viscous flux vector S at the advanced n+1 time step are in general nonlinear functions of the unknown 27 vector q. They are linearlized while maintaining second-order time accuracy as follows. The flux vectors are linearlized by using a Taylor series expansion about q" as follows a E E"n+1= E"+AnAq" + O[(Aq)2] where A = aq F"*1 = F"+BAq" + O[(Aq)2] where B (3.3) aq G"*I = G"+CnAq" + 0[(Aq)2] where C = a aq The above approximation is second-order in time since Aq" = q"*1-q" = O(Ar) and O[(Aq)2] = 2[(Ar)2] The viscous flux vector S is linearized by using the procedure suggested by Steger [291. The coefficient of viscosity y and thermal conductivity k are assumed to be locally independent of Jq ( recall that Jq is the unknown vector in physical domain ), so that elements of S can be expressed in the general form si = J-1, ai where ct is independent of Jq. a Each element is then linearized in the following manner by Taylor series expansion sin+I = s "+J-1{C : )"]}JAg * + 0[(Ar)2 a aJqj Written in vector form, the viscous flux vector becomes S"*I = S"+J-'M"JAq" + O[(Ar)2] (3.4) 28 The matrices A, B, C and M are so called Jacobian matrices and the elements of these matrices are given in Appendix A. After substituting all the linearized terms, Eqs.(3.3- 3.4) into Eq.(3.2), one has [I+ AT( + + -Re-, a (J-IMJ)"Aqn 2 ac an a a = -Ar( aE + aF + --Re-, )" + O[(Ar)2] (3.5) a8 an a a where I is the identity matrix, and Aq" is the delta form unknown to be solved for. Note also that a first-order time-accurate Euler implicit scheme can be obtained by simply replacing Ar/2 on the left-hand side of Eq.(3.5) by Ar [281. 3.1.2 Approximate Factorization Solving Eq.(3.5) directly is very inefficient since it involves the inversion of a very large system resulted from three dimensions. Beam and Warming applied an approximate factorization procedure to reduce the inversion to three onedimensional inversions. They rewrote Eq.(3.5) in the factored form Ar 8A Ar B Ar a + ]"[+ ]"[+ A -Re-' --(J'1MJ)]"Aq" 2 a 2 an 2 a a = -Ar( aE + aF + -G-Re-, as )" + O[(Ar)2] (3.6) a an a a Since the error introduced by this factorization procedure has the leading thirdorder term Ar3 aA aB aB aC +C 8A ) aq" 4 a8 an an a a a ar the second-order accuracy of the difference approximation is not affected. 29 3.1.3 Spatial differencing The spatial derivatives appearing on the left-hand side of Eq.(3.6) are approximated by second-order finite-difference formula. e.g. - 6f + O[(A )2I _ 1-(f-fi.1) + O[(A )2 hence the resulting scheme possesses a block-tridiagonal structure that can be solved efficiently. On the right-hand side, the same second-order differencing is employed for the viscous term, while a fourth-order finite differencing approximation, e.g. af 6 - + O[(AI)4] (fi+2+8fjar-8fj-1+fi-2) + O[AO)41 at12Aq~ is used for the convection terms to keep the convection truncation error from exceeding the magnitude of that of the viscous term. After applying spatial differencing into Eq.(3.6), the spatial derivatives are replaced by difference operators, and the finite-difference form of Eq.(2.21) is then [I+ "[I+B "+ (6,C-Re-'s,(J-'MJ))]"Aq" = -Ar(6,E+6 1F+6cG - Re-16,S)" (3.7) This is the basic Beam-Warming algorithm. The implicit operator has been decoupled to produce an ADI-like scheme, and due to the second-order central-difference operators employed, the algorithm produces block tridiagonal systems for each spatial direction. 30 3.1.4 Numerical Dissipation The Beam-Warming scheme requires the use of numerical dissipations to damp out the spurious oscillations that occurs when dealing with discontinuities. Several different numerical dissipations have been investigated and discussed in [30]. Those employed in the original ARC3D code were proposed by Beam and Warming [28], and Steger [29] and are described as follows: The second-order implicit dissipation operators De = -e1J1VIVeJ D,, = -eJ-Yv A,,,J (3.8) D. = -e1J-VCASJ are added to the implicit operators on the left-hand side of Eq.(3.7) in the , and directions, respectively, and on the right-hand side the following explicit fourth-order dissipation term is added: DE = -EEJ_1V [AvCA2)+ (VA7)2 + (V Ar)2 ]JAq' (3.9) in which symbols A and v are the forward and backward difference operators, e.g. af - --Af, + O[Af] = (fi-1f.) + O[Af] f - VA + O[Af] =-(f1-f1.1) + O[Af] and eI, EE are constant coefficients used to control the amount of numerical dissipations. It was suggested to set EE=At and eI=2EE so that as the time step increases, the numerical dissipations added relative to the spatial derivatives of convection and diffusion remains constant. 31 After inserting numerical dissipations into Eq.(3.7), the finite-difference representation of Eq.(2.21) becomes [I+ A6A+D]"[I2+ 6,B+D ,]"[I+ Ar (6 C-Re-'6,(J-'MJ))+D ]"Aq" = -Ar(6,E+6,,F+6,G - Re-'sS)"+DE (3.10) 3.1.5 Solution procedure For simplicity, the following notation is introduced Le = [I+ 6,A+D] L,7 = [+ .6,,B+D,]"l Lr = [I+ A(6C-Re-'6 (J-'MJ))+Dc]"Aq" R" = -Ar(6 E+6 ,,F+ 6 G-Re-16,S)"+DE where the L's are operators, then the Beam-Warming scheme, Eq.(3.10), can be written as L L,,LAq" = R' (3.11) If one further defines intermediate solutions as Aq' = L rAq" Aq* = L,,Aq* the unknown vector q at the n +1 time level can be computed as follows. LcAq** = R" L,,Aq* = Aq* LcAq" = Aq* (3.12) 32 qn+ = q"+Aq" The Beam-Warming scheme has been shown to be accurate and efficient for simple flow problems (10], however, for more complex flow cases it can be very sensitive to grid resolutions [25]. 3.2 TVD Scheme It is known that finite-difference schemes may produce oscillations when applied across a discontinuity. For instance, the second-order Lax-Wendroff scheme un+ I= iuj -v( u",.-u,' )-Hv(1-v)( u"+-2u"+uj' ), v = aAt/Ax (3.13) for a linear hyperbolic equation au au -- + a-- = 0, a > 0. at ax can produce spurious oscillations near the discontinuity. The Total Variation Diminishing (TVD) scheme is one of the methods designed to eliminate those undesirable oscillations. The notion of TVD was first introduced by Harten [31]. He derived a set of sufficient conditions from which the constructed schemes will not generate oscillations across discontinuities. Since the theoretical background of TVD schemes is rather complex and involved, only relevant basic definitions are given here. The construction of a symmetric TVD scheme and its application to thin-layer Navier-Stokes equations is discussed in the subsequent sections. 33 A numerical scheme, either explicit or implicit, can be written in the following operator form [32] L U"*1 = R u" (3.14) where u"+1 is the mesh function at the advanced time step to be solved, u" is the known mesh function at current time step, and L and R are some finitedifference operators designed for the particular scheme, e.g., for the second-order Lax-Wendroff scheme given in Eq.(3.13) L 1 R = 1 R =1 -vAj+4-Mu2(1-v,)Aj.,A Aj Here and in the following, ey stands for the central difference operator, i.e., Aiu = uj+i - uj. The total variation of a mesh function u' is defined to be TV(u") = .z I u'y1-u"3I = . I A+ u" (3.15) J=- _Wi =Here, the notation - is used to emphasize that the summation is taken over the entire domain. If a numerical scheme applies to a mesh function u, the total variation of the mesh function is non-increasing as a function of time, i.e. TV ( u"+1 ) s TV ( u" ) for all n the scheme is said to be TVD. In practice, the following sufficient conditions derived by Harten [31] are served for constructing TVD schemes: TV ( R u" ) TV(u ) TV ( L u"*' ) TV ( u"') (3.16) That is, if the numerical scheme shown in Eq.(3.13) satisfies the sufficient conditions, Eq.(3.16), the scheme is TVD. 34 It has been proved that a TVD scheme is monotonicity preserving [32], that is, for any monotone mesh function u that increasing (or decreasing) monotonously with time, w = Qu is also monotone, where Q is a finite-difference operator. Since a monotonicity preserving scheme will not generate spurious oscillations, and neither will the TVD scheme. In what follows, the construction of the symmetric TVD scheme for a onedimension hyperbolic equation will be described first, then the scheme is applied to the Euler equations, and finally it is extended to the thin-layer Navier-Stokes equations. 3.2.1 Symmetric TVD Scheme The symmetric TVD scheme was developed by Yee [27] and it is a generalization of Roe [33], Davis [34] and Sweby's [35] TVD Iax-Wendroff scheme. To illustrate the construction of this scheme, a one-dimension linear hyperbolic equation of the form 8u au -- + a - 0 (3.17) at ax is considered. Sweby has shown that the second-order Lax-Wendroff scheme for this equation can be constructed by adding an anti-diffusive flux to the first-order upwind scheme [35], i.e., uji= u1 - vIgAu" - M.(l-v)A.)Aj+%u", if a > 0. i4 + = i - Lj+%u" - v(1+v)j+%A. u" , if a < 0. (3.18) 35 where v = aAt/Ax is the Courant number. Then he introduced a limiting function p and rewrote the above scheme as 11+1 = uj - V 1,+(1-v)[p(rj')/r,*-p(r'j.1)]}ag.gu" if a > 0 (3.19) 114+ = u1 - v{1-%(1+v)[p(rg-,)-p(r-j)/r-j]}Aj,4u" if a < 0 where p is chosen to be a function of the ratio of consecutive cell gradients rj; p, = p(rj), r,+ = Aju/A,.+,u, r; = A+,u/A-u, with the restrictions (sp/rs2, 0sp:2 to satisfy the TVD constraints, Eq.(3.16). The Sweby's scheme, Eq.(3.19), was reformulated by Davis [34] as a LaxWendroff scheme plus an upwind weighted numerical dissipation term, i.e. 14+1 =j - 1/(1+v)A u"- /hv(l-v)Aj u" + [K (rj+)+K-(r- 1)]Aj.,u" - [K (r*j-,)+K-(r-)]Aj.u" (3.20) where K v(1-v)[1-p(rj+)] if a > 0 1.0 ifas0 K-(r) 0 if a <: 0 (3.20a) v(l+v)[p(r- )-1] if a < 0 To avoid determining the upstream direction, Davis then modified the dissipation term as follows. K+(rj+) = jv|(1-|vj)[1-p(r')]/2 K-(rj ) = Jv|(1-|vI)[1-p(r;)]/2 (3.20b) Shortly after that, Roe [331 further generalized Davis' scheme, Eq.(3.20), and suggested a new limiting function Q. Roe's scheme takes the form 36 uj+7= uI - 2v(l+L)A.4u" - %v(1-V)AjU" - Vvj(1-|Vj)(1-Qj.4)Lj_%U" + V j (1-jvj)(1-Qj %)Aj 4U" (3.21) where the limiter Q depends on three consecutive gradients given by QJ+% = Q(Aj.-u/Aj+%u , Ay1,u/Aj+,u) = Q(r*j , r- 1) (3.21a) and 0 < Q < 2/(1- v) 0 < Qyg/r-,i < 2/Ivl (3.21b) 0 < Q3+,/r*j < 2/1&,l to satisfy TVD sufficient conditions, Eq.(3.16). There are a variety of limiter Q satisfing the above restrictions, such as Q(r',r+) = max[O, min(2r-,1), min(r-,2)] + max[O, min(2r*,1), min(r*,2)] + 1 (3.21c) Roe's scheme was later modified and extended to a one-parameter family TVD scheme by Yee [27]. She first rewrote Roe'scheme, Eq.(3.21), into the form uj*1+ Ae(h,'%-hj'g1) = uj" + A(1-e)(hj,-hj",) (3.22) hj+ = [a(up+1+uj)-{Aa2Qj+&+IaI(1-Qj+4)}Aj+4u]/2, A = At/Ax (3.23) and then simplified the numerical flux to hj+4 = [a(uji+uj)-ala(l-Qj, .+ u]/2 (3.24) Note that the term Aa2Qj+%Aj+%u/2 was taken out from Eq.(3.23), but it makes no difference as long as the limiter function Q is chosen such that the scheme satisfies the TVD sufficient conditions. 37 For nonlinear hyperbolic conservation law of the form + af(u) = 0 (3.25) a t a x the numerical flux, Eq.(3.24), was modified to h = [(f ,i+fj)-w(a , )(1-Qjg)Aj 4u]/2 (3.26) where f=f(uj), A + u=uj,-u+ , a(u)=df(u)/du is characteristic speed and aj+ { (fj+rffj)/j+%u , - 0 (3.26a) a(uj) , Ay+u 0 The function F(a+,) that is sometimes referred to as the coefficient of numerical viscosity, is defined as { |zI , IzI 3 2 b 'i2(z) = {(2+z2)/2c ,z< (3.26b) to avoid the entropy violating problem [31], where E is a parameter that is described in Ref.[36]. In order for the above scheme, Eqs.(3.22, 3.26), to be TVD, that is, satisfying the TVD constraints, Eq.(3.16), the following conditions obtained by Yee [37], must be satisfied 0 < Qi+ < 2 0 < Qj+ /r*j < 2/[v(1-9)a,.g,] - 2 0 < Qi+4/r-;+1 < 2/[L.(l-e)Ia,+1I1 - 2 (3.26c) v Iaj+, I< 2/[3(1-e)] If one selects the following limiter that satisfies the above constraints Qi% = minmod(Aj.u , Aj ,u) + minmod(Aj ,u , Aj 1,u) - Au+u (3.27) then from Eqs.(3.26-3.27) Yee's numerical flux can be expressed as hj+%= [fj + fj + j+gJ/2 (3.28) O+ = *(a,+)(gj+gj+1) - 2,&4.,u (3.28a) 38 where the minmod function minmod(xy) is defined as 0 , if x y have opposite sign minmod(x,y)= 9 and gj = minmod(t-.,.u , j.,u) (3.28b) The above scheme, Eqs.(3.22, 3.28), form a class of TVD schemes with the numerical dissipation terms centered and hence are often referred to as symmetric TVD schemes. It is noted here that the scheme is TVD for onedimensional nonlinear scaler hyperbolic conservation laws, and is also TVD for constant coefficient hyperbolic systems. However, in applications, it is formally extended to more complex equations, such as the Euler equations, the thin-layer Navier-Stokes equations, and evaluated by numerical experiments. 3.2.3 Extension to Euler equations Extension of the scaler TVD scheme to systems of conservation laws is not unique. Here, one follows the method proposed by Yee [38,39]. The method is accomplished by defining at each point a "local" system of characteristic fields and then applying the scheme to each of the scalar characteristic equations and hence is referred to as the "local characteristic approach". Consider the inviscid part of the thin-layer Navier-Stokes equations, Eq.(2.21), aq aE aF aG -- + - + -- +--- = 0 (3.29) at a an; a 39 Recall that the Jacobian matrices of E, F, G are A, B, and C, respectively. Denote the eigenvalues of A, B, and C as (al a ai '), a~ia~ia) (4i,agaaa") and Rt, R, R, as the eigen-matrices whose columns are eigenvectors of A, B, C, respectively, and R-'E, R',, R-, as the inverses of the corresponding eigenmatrices. The explicit form of these eigenvalues and matrices are given in Appendix B. In the -direction, let q,,, be some symmetric average of %, and q,, e.g. qj+% -= 0.5(q 4+q,,W) or Roe's average % + = (P gjq*+ p 4j1,jqy 1W)/(PWh + Phy ) (3.30) and denote a',+, R,, R1j+% as the quantities a',, R,, R71, evaluated at qAk. The difference of the characteristic variables in the local -direction are defined as aj+4 = R+'j g(J ,,Wqj +,-Jgq+)/[(Jj, +J J)/2] (3.31) One can also define the corresponding parameters for the r7 and directions in a similar manner. Using these notation and applying Eq.(3.22) to each equations, a oneparameter family of TVD scheme for the system of equations can be written as q + Ae[(Ej";* j-E _* F + " F") +(Gj" -Gj')] = qJk - \(1-9)[(Ej",sW-E _4,)+(F "k -Fjl_%, +(GjL+ 4-Gj"0.%)] (3.32) 40 where A=Ar since A =Ar =A =1 and the numerical flux functions obtained from Eq.(3.28) are given by Ej+* = [Ej j+Ej+Ryg4Ds]/2 Fj,k+, = [Fgky+Fj,k+lI+Rk+ 4Dk+ ,]/2 (3.33) Gj4.,.% = [Gjjo+Gj,+RaseggD,]/2 the expression of the elements of Dj.,. are obtained from Eqs.(3.28a-3.28b) and are given as 0 += t(aj+,)(egj, -2aj+,) (3.33a) g'. = minmod(4x. ,a.+%) (3.33b) Likewise, one can define Dk+4, II+4* 3.3.6 Linearization and ADI decomposition In order to solve Eq.(3.32) efficiently, the left-hand side of Eq.(3.32) is first linearized by the procedure described in Section 3.1.1. The resulting scheme written in delta form is [I+ Ae[(Hj+ %*+-Hjg +(Hj+j-Hjk-.,)+(Hjg g-H p.g)]Aq" = -A[(Ej +,4-Ej.-,1) + (Fj+,.j-Fk,1) + (Gjg,+4-Gj%)]" (3.34) where I is a 5 X 5 identity matrix and the H operators are given by: Hj+%,, = [A +n s,]/2 H.,k,, = [Bjk+1 j,k+4,I]/2 (3.34a) Hjgj+ = [C +~ogj,I+]/2 in which nj+.k = {R diag[p'-2,y(a')] R-1}+qA+% 41 Qk+m, = {R diag[p'-2k(a')] R-1}k+4Ak+4 (3.34b) aj4, = {R diag[p'-2,Y(a')] R-1,,.A and flj+ = *(aj 4)(gi +gjj )/ct'j , += k k+-Jdk+ek+0/aik+A (3.34c) Si.. =+ *(a',.j,g(g e,ig1a'j , Here, the expression diag(z) denotes a diagonal matrix with diagonal elements z'. To calculate the o's at every time step is quite costly. For steady-state applications, the temporal accuracy is not important. One can use a spatial firstorder implicit operator (that is, by setting the limiters equal zero); i.e., by redefining the n's as nj+%,v = {R diag[-*(a)] R-1})+,+ ak+%4 = {R diag[-*(a')] R-1}k+%Ak+% (3.35) njW+4 = {R diag[-i(a')] R-1}1,A1+ Since the temporal accuracy is not important, the computation can be further reduced by simplifying the n's diagonal matrices, i.e. n+v.*ji = {diag[-maxt(a')]},A+, ok+v4 = {diag[-max(a)]}k+%Ak+% (3.36) nj., = {diag[-maxY(a')]},Asa If one applies the approximate factorization procedure discussed in Section 3.1.2 to Eq.(3.34), the ADI form is obtained, i.e. [I + Ae(Hj+ %,-H.%,k)][I + ) (Hjk.,.-Hjk-&,)][I+ A(H ,.%-Hj.%)]Aqn = -A[(Ej+%*k-Ej_. ) + (Fjk-+%I-Fj,k-,) + (Gj,+ -G .%)]" (3.37) 42 With the use of Eq.(3.33) and Eq.(3.34a), the above equation can be expressed in delta form as follows [I+ Ae6,A+D,]"[I+ Ae6,,B+D,,]"[I+ As6,C+D,]"Aq" = -A(6fE+6,,F+6G)"+DE (3.38) where D= D,= (3.38a) D= DE = -,+R-R.i+RR + s-14%4-R1_%oj.%) (3.38c) 3.3.7 Application to thin-layer Navier-Stokes equations To apply the TVD scheme to the thin-layer Navier-Stokes equations, Eq.(2.21), is similar to that to the Euler equations [38-39]. The viscous term that not appears in the Euler equations is first central differenced and linearlized as discussed in Section 3.1.1 and then added to Eq.(3.38). The resulting finitedifference equation is [I+ Ae6,A+D,]n[I+ +e6,,B+D,,]" xe(6,C-Re-'6,(J-'MJ))+D,]"aq" = -A(6fE+6,,F+s G - Re-s 6S)"+DE (3.39) It is noted here that A=At since A =Ar =A =1, and e is a parameter in the range of zero to one; Euler explicit or implicit scheme is obtained by setting 9 = 0 or 1; if e =1/2, it is a second-order implicit scheme. However, the order of accuracy is degraded due to the simplifications that have been made in the derivation of the scheme. It is also noted that if 9 =1/2, Eq.(3.39) is similar to 43 the Beam-Warming scheme, Eq.(3.10), except that the numerical dissipation terms are replaced by the ones derived from TVD constraints. In fact, the TVD scheme can be viewed as a central difference scheme with a "smart" numerical dissipation which is an automatic feed back mechanism used to control the amount of numerical dissipation [38]. 3.3 Flow Solver Applications of the original version of ARC3D and its axisymmetric version [40] for the simulation of transonic turbulent flow past a projectile model showed that the codes implemented with the Baldwin-Lomax turbulence model can give surface pressure distribution in excellent agreement with measured data if good grid networks are provided to the Navier-Stokes codes for the simulation [25]. Numerical experiments also showed that the Beam-Warming solution algorithm is very sensitive to grid characteristics in complex flow regions. In order to improve the robustness of the codes, a symmetric TVD scheme was first extended for modifying the axisymmetric code and investigated for the projectile aerodynamics simulation; as reported in [25,41], the modified code is indeed less sensitive to the grid characteristics in complex flow regions. Later on, the original ARC3D was modified to implement the symmetric TVD scheme. As have been mentioned before, from the viewpoint of numerical implementation, the only difference between the Beam-Warming scheme and the TVD scheme are the numerical dissipations. Hence, to implement the TVD scheme into the original ARC3D, one simply replaces the dissipation operator, Eq.(3.8), and the 44 dissipation terms, Eq.(3.9), with those designed for the TVD scheme, i.e., Eqs.(3.38a, 3.38b). However, the latter is more complicated than the former. New subroutines for computing the TVD dissipations had been written and the original ARC3D was then modified into a TVD thin-layer Navier-Stokes code [26]. This TVD code has been applied to the projectile flow problem, and the results showed that the code was indeed accurate and robust. Hence, it is adopted in this study for developing the multi-block computational method. CHAPTER IV THE FLOW PROBLEM AND SOLUTION STRATEGY A transonic turbulent flow of Mach 0.8 past a wing-fuselage configuration shown in Figure 4.1, is adopted as the model problem to develop the multi-block computational method. The flow problem considered is indeed a complex one, involving many viscous phenomena, such as viscous sublayer, shock boundary layer interaction, flow separation, as well as complicated boundary geometry. There are wing surface pressure measurements available for assessing the accuracy of the numerical solutions. The flow field is assumed to be symmetric with respect to the symmetry plane of the wing-fuselage configuration and hence only half of the flow domain is required for the simulation. 4.1 The Flow Problem and Preliminary Work The wing-fuselage configuration shown in Figure 4.1 is a model of a nextgeneration 150-passenger transport. The wing, optimized for Mach 0.77, has an aspect ratio of 5.17 based on the averaged chord length, wing leading edge sweep of 30*, and an advanced supercritical airfoil. This wing-fuselage configuration is used as a model to develop the computational methods for propfan powered aircraft [42]. A rather coarse one block H-grid for the configuration was provided by the NASA Ames Research Center. The flow field is assumed symmetric with respect to the z = 0 plane and hence only the z > 0 half-space is 45 46 Ilk ~ i': Figure 4.1 Wing-fuselage configuration .,wi% 47 considered for the simulation. Figure 4.2 shows an overview of the physical domain and its computational domain. In the following, the integers j, k and 1 are used to denote the grid points in the , 7, and directions, respectively. Here, is along the streamwise direction, t is along the wing spanwise direction, and is about normal to the solid surface. The dimensions of this original grid network are jmax = 101, kmax = 45, and lmax = 29. The solid surfaces are transformed into part of the 1=1 plane with the shape H, as shown in Figure 4.3; k =1 and k = kmax are symmetric planes; while j =1, j =jmax, and I= lmax are artificial outer boundaries. To further illustrate the configuration quantitatively, one takes the maximum diameter of the fuselage, D, as a reference length. The fuselage is about 7D long; the wing span is about 3.3D; the chord length at the wing root is about 1.4D; the diameter of the sting is about 0.075D; the distance of the outer boundaries far upstream and downstream are about 15D, while that of the lateral side is about 10D. These relative lengths and the j, k, I triplet in physical space are shown in Figure 4.4. Note also that the first spacing next to the solid surface is in the range of 0.008D to 0.026D. The 101x45x29 H-grid one-block system was first used for a preliminary study. Because of the specific grid topology, the thin-layer Navier-Stokes code simply can not be applied directly without modification for turbulent flow simulation. As can be seen from Figure 4.3, only part of the 1=1 surface is the solid surface. Extensive modifications were made for proper Jacobians and metrics calculations, boundary conditions, and eddy viscosity computations. Also, 30 grid points that extremely close to the solid surface were added to the original 48 9-2C cF E L JK45 z- J-1 D A 00 Qp Y GH . C Figure 4.2 An overview of the original grid 49 (a) Physical space 0 U v PH NOSE FUSELAGE TAIL STING 0 R F WING K L WING E 0 R F NOSE FUSELAGE TAIL STING 0 S (b) Computational space Figure 4.3 The 1=1 surface of the original grid K L F E 0 z S, T w 8D I 1,23.,1) 3.GD ..C(1.,45.,29) 37 ,23,1) E-15E K L 76 23., ,29) 0 p J -t,2131, 1 (13,l.,l) 7D Q CA(90 1 )10D Z37D' 15D F D GH (101, 45,29) y B 1OD 0.075D x (101,1l,29) 51 grid to resolve the thin viscous-layer; and two mirror image surfaces were also added for ease in symmetric boundary condition treatment. The 101x47x59 Hgrid network was then provided to the modified thin-layer Navier-Stokes code for the flow simulation. As shown in Figure 4.5, the computed results are in poor agreement with experimental data. This seems to indicate that the number of grid points on the wing surface, 18 in the spanwise direction and 40 in the chordwise direction, is not fine enough to resolve the complex flow phenomnia. Also, modifications on the thin-layer Navier-Stokes code may be needed before the one-block H-grid topology can be used for this complex turbulent flow simulation. A flexible and efficient multi-block technique is more desirable to numerically simulate a complex geometry turbulent flow problem, such as the wing-fuselage aerodynamics. It is also noted here that the calculations of this one-block approach were performed with the Cray X-MP/48 of the NASA Ames Research Center. 4.2 Solution Strategy For flow past a complex configuration, such as the wing-fuselage combination, to generate an acceptable single grid for the entire flow field is rather difficult. It is much easier to divide the flow field into several subregions and then generate a good grid for each subregion. Hence, it is natural to develop a solution method based on the multi-block grid network by treating each block as an independent flow problem. Each block is then solved by the same thinlayer Navier-Stokes code with the information between blocks transferred via the interface boundary conditions. 52 /P 0 M - M - 0.8 -000 00 o0 0 o - 0 a.-.4 uym 6- I ( / / / / I / / / / / / / -I /1 Figure 4.5 Computed results of the one-block approach CP. I 53 For the wing-fuselage flow problem, an effective and yet simple block-byblock computational method is proposed. In this method, the entire flow field is divided into six contiguous blocks and each block is partially bounded by a solid surface. Each block is then treated as an independent flow problem, which is solved by the same thin-layer Navier-Stokes code, with the interface boundary conditions updated at every time step. This method does not need to worry about patching solutions between different equation zone as in the zonal equation approach, and it is easier to deal with the interface boundary conditions than that of the Chimera overset grid approach. However, in the same manner as the zonal methods, poor interfacing can introduce errors at the interface boundaries and consequently propagated into the solution field. The error introduced can depend on the solution behavior and the grid characteristics around the interface as well as on the interpolation technique employed for updating the interface boundary conditions. For a successful multi-block computational method, it is then important and critical to minimize the associated interface error and its effect on the solution accuracy. To attack the wing-fuselage flow problem by using the proposed method, a flow simulation package has been developed which includes four independent computer codes; a surface grid generation (SG3D); a volume grid generation (GG3D); a multi-block flow solver (WF6B); and a plotting program (PP3D). They will be briefly discussed in the subsequent chapter; however, a detailed description can be found in a reference manual [43]. Most of the figures presented in this dissertation were produced by the plotting program PP3D. 54 4.3 Six-Block Grid Topology For a multi-block computational method, the definition of a good composite grid topology for a complex configuration turbulent flow simulation is more involved and difficult than that for an inviscid flow or a laminar viscous flow simulation. For instance, the solid surface has to be transformed onto an entire plane in computational space so that the thin-layer approximation can be applied easily; the Baldwin-Lomax turbulence model requires the normal distance of each grid point to the solid surface; and the grids should be able to cluster points close to the solid surface to resolve the thin viscous-layer. Hence, for ease in eddy viscosity calculation and in applying the thin-layer approximation, a special 6-block grid topology was defined and investigated for the wing-fuselage configuration. The main idea of dividing the flow field is to isolate the wing by a cone-like boundary surface that sits between the wing and the fuselage, as shown in Figure 4.6. The isolated wing was divided into 2 blocks and the rest of the flow field was divided into 4 other blocks. By doing so, the flow field is divided into six contiguous blocks such that each block is partially bounded by a solid surface, and the solid surface of each block is transformed onto = 1 plane in the computational space. An overview of the topology is shown in Figure 4.7 while Figure 4.8 gives side and cross-sectional views. The physical and computational space of each block is shown in Figure 4.9. The selected grid topology has the clear advantage that each block can be solved by the same thin-layer Navier-Stokes code; however, the choice has its deficiences in the block grid generations. For instance, for the wing blocks, Block 55 zY dividing surface z Y ividing surface Figure 4.6 The dividing surface isolating the wing 56 -o I NAV CC / c0 C-C UtJ Figure 4.7 Six-block grid topology (overview) 57 Figure 4.8 Side and cross-sectional views H S T N F BOC 1\C /BLOCK 5 x A B K Z G M L x U z H,S,T,N b' BLOCK 3 BLOCK 4, R,Q BLOCK 2 BLOCK 5 b 0 G, m B,K D,P 0,1 58 H E A B GB BLOCK I Figure 4.9 The transformation of the six block grids z y F G Ht C AE DD B C D (a) Block 1 H z y N G N CI B b CL 2 -L K (b) Block 2 G - - - - - - - -B- M BLOCK 59 S N H~f )I/ fitI \_e QR C C L bL (c) Block 3 Y H N T R C L BLOCK 4 (d) Block 4 Figure 4.9 ( continued ) B BLOCK 3 s T N C L AH s T (e) Block 5 z y -- x x z U y K L K (f) Block 6 Figure 4.9 ( continued ) 60 0 z N U P K BLOCK 6 C L p-- - -- -8125 H 00 N C L 61 3 and Block 4, shown in Figure 4.9(c,d), four of the six bounding surfaces in the computational domain, CLRQ, LNTR, RTSQ, and QSHC form a nearly flat surface in the physical space. Therefore, special measure and care in grid generation are required to overcome the difficulty resulting from the drastic change in the domain transformation. Also, for the nose block, Block 1, shown in Figure 4.9(a), the boundary plane AEJF in the computational domain degenerated into a single line in physical space and care must be taken when dealing with the boundary conditions. CHAPTER V GRID GENERATION For the wing-fuselage flow problem, the flow domain has been divided into six blocks, however, it is still a great challenge to generate acceptable grids for each block. To generate grid networks for the six-block wing-fuselage flow problem, surface grids must be obtained first. A special surface grid generation code (SG3D) was developed to utilize the original H-grid network for constructing the required surface grids. The resulting surface grids will be presented in section 5.1. Generating volume grids for the six-block flow domain, it is required to develop an adequate grid generation code. Various grid generation methods [44-51] have been studied in that regards. The main concern is that the code can produce grid networks that satisfy the required grid characteristics, such as: (a) grid lines in different directions must be nearly orthogonal to each other; (b) grid lines in c-direction must be about normal to the solid surface; (c) = constant grid surfaces must be clustered near the solid surface; (d) grid lines must be continuous across the interfaces; and (e) the blocks must be contiguous (no holes or overlapping). A combination of two-boundary algebraic and elliptic grid generation methods will satisfy these requirements best. Accordingly, a general purpose 3-D Grid Generation code (GG3D) that based on these two methods was developed and applied to generate the 3-D grid networks. 62 63 A detailed description of GG3D can be found in [43]. The techniques employed to generate the six block grids are described in Section 5.2 and the resulting grid networks are presented in Section 5.3. 5.1 Surface Grid Generation To generate grid networks for the six block grids, surface grids must be obtained first. A two-boundary algebraic grid generation method is used to generate volume grids; consequently, only the solid surface grid (1=1) and the cooresponding outer surface grid (1=imax) are required. The other boundary surfaces are then defined by the straight lines that connect the boundary curves of the solid and outer surface grids. A special surface grid generation code (SG3D) has been developed that uses the original H-grid to construct each surface grids by cubic splining, transfinite interpolation, strentching, shifting and dividing. The resulting surface grids obtained are shown in Figure 5.1; the dimension of these surface grids are 28x23, 40x12, 40x35, 40x35, 40x12 and 26x23 for Block 1-6, respectively. It is noted here that the outer surfaces are defined according to the grid topology and will have strong effect on the grid characteristics of the block grid to be generated. Hence, the boundary curves of each outer surface should be defined carefully. For block 1, a cone-like shape that is resemble to the solid surface is desired for the outer surface. As shown in the grid topology, Figure 4.7, the boundary curve GHI is defined such that the bounding surface lies in the middle of the wing leading edge and the fore-body of the fuselage. Similarly, the outer surface of block 6 is defined such that bounding surface of the boundary 64 NA (a) Solid surfaces Figure 5.1 The 6-block surface grids lapr 6 65 fir (b) Outer surfaces Figure 5.1 (continued) N 66 curve MNO lies in the middle of the wing trailing edge and the aft-body. The boundary curves Hb'N of block 2 and Ha'N of block 5 are also defined such that the boundary surfaces Hb'NLC and Ha'NLC lie in the middle of the wing and the fuselage. For the outer surfaces of block 3 and block 4, the boundary curve HN, that is defined as a straight line, is located such that the interface surface of block 3 and block 4, i.e., surface HRCQLN, is not incline to either upper side or lower side of the wing. In other words, the plane spanned by HCLN should pass throught the wing tip, and hence divides the wing into upper wing and lower wing. If this were not the case, one would have difficulty to generate an usable volume grid for either block 3 or block 4. 5.2 Volume Grid Generation 5.2.1 Two-Boundary Grid Generation The two-boundary grid generation technique used is based on Hermite transfinite interpolation [46-48]. Two surface grids are specified and the interior grid lines are connected by curves such that boundary orthogonality is enforced. The other four boundary surfaces are then defined using straight lines. The distribution of interior surfaces is defined by a suitable clustering function such that more grid points are clustered near the solid surface as it is usually required. A description of clustering functions is given in Appendix C. Denoting the position vectors of the two specified boundary surfaces as rQ,,0) = r1(C,q) and r(Q,r,1) = r2(Q,q) (5.1) and the derivatives with respective to as 67 DR, = ar(Q, i,0) and DR2 = ar(Q,,1) (5.2) a a the following functional form, based on Hermite transfinite interpolation, may be written for the connecting curves: r(,v, ) = r(4,t)hj( ) + r2(,i)h2() + DRih3(g) +DR2h4( ) (5.3) To satisfy Eq.(5.1) and Eq.(5.2), the following must be true: hl(0) = 1, hI(1) = 0, ah1(0) = 0, ahb(1) = 0 a a h2(0) = 0, h2(1) = 1, ah2(0) = 0, ah2(1) = 0 h3(0) = 0, h3(1) = 0, ah3(0) = 1 ah3(l) = 0 (5.4) h4(0) = 0, h4(1) = 0, ah4(0) = 0, ah4(1) a . a . Hence, a cubic polynomial in for hj( ) satisfying conditions given in Eq.(5.4) can be determined. They are h= 2 3-3 2+1 h2 =-2 3+3 2 h3 = c3-2 2+ (5.5) h4 = _-2 Here, is a normalized point distribution, and is usually defined by a clustering function as shown in Appendix C. The expressions for the derivatives, DR, and DR2, appearing in Eq.(5.3) are the key elements used to enforce the boundary orthogonality and are derived as follows. 68 A vector normal to the surface r = 0 that points into the spatial domain is given by n - ar1(,y) ar1(Q,7) a a 1 where x denotes the vector cross product. A vector that tangents to the connecting curve at r =0 can be defined as e = ar(Q,2,O) a. Hence, in order for the connecting curves to be normal to the r = 0 surface, the cross product of the vector tangent to the connecting curves e and the surface normal n must be zero, i.e., e x n = 0 at = 0 (5.5) or written in component form ay ax ay ax ay az az ax az ax a a 1 a a a 1 a a n a a a n + ax ay ax ay ax az ay az ay az a a? ae a a1 a an ae ae a n + ax ax az ax az ay az 8) 8z ay )]k a a7 ae ae a17 a( a7 ae ae an =0 (5.5a) From the above equation, it can be seen that the connecting curves will be perpendicular to the surface r = 0 if the following conditions are satisfied at r = 0: DX= - ( ay az ay az a at ac ae an 69 DY, - ay = CK,(,o)( ax az ax az (5.7) a at a a an az ax ay ax ay DZ1 = = CK1(,n)( an a a an Here, CK, are arbitrary functions introduced to control the shape of the connecting curves, and if set to zero the curves become straight lines. In a similar manner, one can show that the connecting curves will intersect the - = 1 surface perpendicularly if DX2 = x -CK2.Q,7)( ay a ya an a a an DY = = CK2(ci) --ax az -z) (5.8) aa ay ax ay DZ2 _ az = -CK2(, x)( an a a an are satisfied at g = 1. By substituting Eqs.(5.7-5.8) into Eq.(5.3), one obtains the following expressions for the resulting volume grid with boundary orthogonality enforced: x(Q,n,,) = x( ,ii)hl( )+x2(,7)h2( )+DXlh3( )+DX2h4( ) y(,n, ) = yl( ,7)hl( )+y2( ,7)h2( )+DYlh3( )+DY2h4( ) (5.9) z(,q, ) = z,( ,t7)h,( )+ z2(,tI)h2( )+DZh3( )+DZ2h4( ) It is noted here that the derivatives on the right hand side of Eqs.(5.7-5.8) are discretized by proper finite- difference fomula, e.g., central difference ay1 _ yl(k+ 1)-y(k-1) + O(An2) an n(k+ 1)-a(k-1) and e and q are normalized as follows: 70 0(j) = (j-1)/(jmax-1), j=1,2,...,jmax n(k)= (k-1)/(kmax-1), k=1,2,- - -,kmax (5.10) The CK's appearing in Eqs.(5.7-5.8) have a strong effect on the shape of the connecting curves. In fact, they control the depth of the perpendicular part of the curves. Values for the CK's are usually specified as constants by trial and error. A shape function was devised to modify the CK's so that different values can be specified at different j, k locations. The shape function is of the form CK = CKI*cosw{CPJ[(j-1)/(jmax-1)]+CSJ} *cosr{CPK[(k-)/(kmax-1)] + CSK} (5.11) where CKI is a constant to be specified, CPJ and CPK are used to control the frequency while CSJ and CSK are used to control the phase of the relative cosine function. In doing so, values for the CK's become functions of j and k, and can be smoothly distributed on the surface. For example, if one specifies CPJ=0.5, CSJ = 0, then CK = CKI when j =1 and CK decreases as j increases; or if one specifies CPJ =1, CSJ = -0.5 then CK increases from zero to CKI and then decreases to zero as j varies from 1 to jmax. Figure 5.2 shows some examples to illustrate this shape function. 5.2.2 Elliptic Grid Generation The technique and subroutines employed are based on the elliptic grid generation system of Thompson [49-51]. The volume grid is constructed by forcing the Cartesian coordinates of interior points to satisfy specially devised Poisson equations. In the Poisson equations, three control functions for each curvilinear direction are set up to control the grid point distributions. The system 71 CPJ-0. 5, CSJ-0 ,CPK-0,CSK-8 CPJ-0,CSJ-0,CPK-.6,CSK--.4 CPIJ-.8, CSJ--O. 4, CPK-.5,CSK--.5 CPJ-0, CSJ-e, CPK-8,CSK-0 CPJ-.5,CSJ--.5,CPK-.5,CSK--.5 CPIJ-1,CSJ--.5,CPK-.5,CSK--.5 CPJ-.6,CSJ--.3,CPK-.6,CSK--.3 CPJ-1,CSJ--.5,CPK-1,CSK--.5 Figure 5.2 Shape function for CK's 72 of Poisson equations are solved by point SOR (success overrelaxation) iterations. The user-provided initial volume grid is used not only to start the SOR iterations but also to set up the control functions. The elliptic grid generation system is based on the Poisson equations of the form V'= (vs.v )Pi V = (vn -vn)P2 (5.13) V2= (v .v)P3 where the P's are so-called control functions which can be specified to control the grid smoothness and/or orthogonality. Since the actual computations are performed in the computational domain where , Y7 and are the independent variables, with r = (x,y,z)T as the dependent variables. The poisson equations, Eq.(5.13), are then transformed to computational domain. The transformed equations are (V. 2 r+(I.V a2 r V. 2r a2r8r a~r (v 7vg) a + (v7.v) + (v vI)- + agoaano agag a2r a2r a2r (v .v-,) + (v7-v-)- + (v .v7) - + aga ana aga (v-v)- + (vn-.vg) + (v-v)- + ar ar ara (v-v )P Ia +(Vn-Vn)P2 a+(VV,)P3 - 0 (5.14) a a?7 a The finite-difference equations of Eq.(5.14) are obtained by replacing the derivatives with central differences and these equations are then solved by point SOR iteration procedure. 73 If the control function P's are evaluated from the initial grid, the three components of the elliptic grid generation system, Eq.(5.14), provide a set of three equations in which the control functions P's are treated as unknowns and the values for the position vectors are obtained from the initial grid. The transformed poisson equations, Eq.(5.14), can be rewritten as (V -V )Pj a + (V17- V1)P2 ar +(v,-V,)P3 arar a1 ar -[( V ) a2r + ( V O a') + (V2.V,) + aga alna asa 8~r 8~r a~r ( -V 7) ar + (vi.vE) + (v .v7) - + ar ar ar (v .v) + (v97.v7) + (v -v,) + (5.15) aa a1W aar which can be solved simultaneously at each point for the three control functions, Pi, i= 1,3. If the control functions are computed solely based on the above equations, and the elliptic system, Eq.(5.14), will reproduce the same grid as the initial grid which is of trivial interest. However, if the computed control functions are smoothed before they are applied into the elliptic system, a smoother grid network can then be produced. A variation to the evaluation of the control functions is obtained by assuming that the grid lines are orthogonal to each other, and Eq.(5.15) is reduced to (v -v )P - ar +(v"-v)P2 ar +(vr.vO)P3 ra a17 as ar ar a~r -[(v -v) + (v.v) + (v .v ) -] (5.16) aa8 a017 asa 74 The control functions evaluated from these equations will force the elliptic system to produce a more orthogonal grid. 5.3 The Six-Block Grid A set of volume grids with 60 grid points in the 1-direction was first generated by two-boundary grid generation with grid lines set normal to the solid surface whenever possible. The first spacing off of the solid surface was set to about 2x10-5D, which is considered fine enough for turbulence simulations. Figure 5.3 shows the 3-D close-up plots of the six block grids with cut-offs to show the details of the interior surfaces. These grids were generated from the solid and outer surfaces described in Section 5.1. The other four boundary surfaces were then defined by straight lines that connect the boundary curves of the solid and outer surfaces. Sixty grid points were distributed on each of the connecting straight lines by using the hyperbolic tangent clustering function with the first spacing of 2x10-5D and the last spacing of six times of the average spacing, and hence, more grid points were clustered near the solid surfaces. To generate grid network using two-boundary grid generation, one has to face a difficulty in specifying a proper value for the parameter CK, Eqs.(5.7-5.8). In general, the value for CK is specified by trial and error, however, its variation on a surface can be justified by the use of the shape function, Eq.(5.11), according to the geometry of the grid to be generated. For the six block grids, the grid orthogonality on the outer surface was not enforced, that is CK2=0, to avoid the possibility of grid overlapping. On the other hand, the grid orthogonality on the 75 1z (a) Block 1 - Figure shows the surfaces 1=1, 1=50, k=5, k=15, j=21 z Y x (b) Block 2 -- Figure shows the surfaces 1=1, 1=50, k=3, k=13, j=17 Figure 5.3 Six block grids -- 3D close-up with cut-off 76 z Lx 1I (c) Block 3 -- Figure shows the surfaces 1=1, 1=50, k=1, k=15, k=30, j=10, j=18, j=32 z x 1 M (d) Block 4 -- Figure shows the surfaces 1=1, 1=50, k=1, k=15, k=30, j=6, j=14, j=32 Figure 5.3 (continued) 77 z Yj (e) Block 5 -- Figure shows the surfaces 1=1, 1=50, k=3, k=13, j=17 Y (f) Block 6 -- Figure shows the surfaces 1=1,1=50, k=5, k=15,j=1,j=10 Figure 5.3 (continued) 78 solid surface is desirable. From numerical experiments conducted the following was one of the best to set up the CK, distribution by the shape function for each block. Block 1 - CPJ=0.5 CSJ=0 CPK=O CSK=O Block 2 - CPJ=1 CSJ=-0.5 CPK=0.5 CSK=O Block 3 - CPJ=0.8 CSJ=-0.4 CPK=1 CSK=-0.5 Block 4 - CPJ=0.8 CSJ=-0.4 CPK=1 CSK=-0.5 Block 5 - CPJ=1 CSJ=-0.5 CPK=0.5 CSK=O Block 6 - CPJ=0.5 CSJ=-0.5 CPK=O CSK=O These distributions are shown in Figure 5.4. For block 1, CK, decreases as j increases and becomes zero when j =jmax. This is reasonable since the j =jmax surface incline about 45 degrees to the solid surface. Similar argument can be addressed for the other blocks. To further improve the smoothness of the grid network, these algebraic grids were then used as initial grids for elliptic grid generation. Smoother grids were obtained for block 1, 2, 5, and 6. However, it failed to generate elliptic grid for either block 3 or 4. This may due to their special geometries such as three boundary surfaces form nearly a flat plane. Figure 5.5 illustrates the difference between algebraic and elliptic grids. Notice that the elliptic grids are smoother and more orthogonal than the algebraic grids. These six block grids were further modified to include mirror image planes for imposing symmetric boundary conditions in the Navier-Stokes simulations. The resulting grid network was then used as a basic grid in this study and is referred to as the "base grid" hereafter. The dimensions of this base grid are 28x25x60, 40x13x60, 40x35x60, 40x35x60, 40x13x60, and 26x25x60 for blocks 1-6, 79 BLOCK I RLIC f3 ~ RLC rWV BLOCK 2 j-jmax kka BLOCK 4 BLOCK 6 Figure 5.4 Shape function for the six-block grid k-kmax J-JBx m nrv 3 I J- jmax k-mkmax k-kmax j- jmOX x k - k max Bi nCK 5 BLOCK 6 6 j-jmax k-kmaxj- jmx 80 (a) Block 1 - k=13 Figure 5.5 Elliptic and algebraic grids z Block 1, k-13, Algebraic Grid z Block 1, k-13, Elliptic Grid 81 z I I- A F/ 7///// / / Block 6, k-13, Algebraic Grid z Block 6, k-13, Elliptic Grid (b) Block 6 -- k=13 Figure 5.5 (continued) Z///, F111 / / / 82 respectively. Figure 5.6 shows some detailed plots to illustrate the grid characteristics. Figures 5.6a-d show that the grid lines are about normal to solid surface whenever possible. Figures 5.6e-f show that some k=constant surfaces are quite skew for block 2, 3, 4, and 5. Figure 5.6g shows the grid distribution at k =1, 12, and 18 wing cross-sections. In order to investigate the effect of grid resolution on the flow solution, a series of grid refinement studies have been performed. Several locally refined grids (near the wing tip or near the solid surface, etc.) were generated and used to study the possible effect due to the grid characteristics. A set of "refined grids" was also generated by doubling the grid points in the wing core-wise direction. The dimensions of these refined grids are 28x25x60, 79x13x60, 79x35x60, 79x35x60, 79x13x60 and 26x25x60. Figure 5.7 shows the solid surfaces of block 2 and 3 of this refined grid. Note that the grid lines were generated by using a cubic spline so that the solid surfaces would not be distorted. 83 z j-1 or 40 Block 1 k-13 Block 3 or 4 k-35 'h ~//A'$"' K (I// <. / / // / Block S k13 (a) Block 1, 3, 6 Figure 5.6 Some detail plots for six-block grids J-40 or 1 (b) Close-up to show the boundary orthogonality Block 3 84 z Black 4 Block 5 J-20 (c) Block 2, 3, 4, 5 Figure 5.6 (continued) B lock 2 i-20 (d) Close-up to show the boundary orthogonality Z 85 Y x k-6 1-1-35 (e) Block 3, 4 -- k surfaces are quite skew k-18,13 1 3 8 1ock 2 I z (f) Block 2, 5 -- k surfaces are quite skew Figure 5.6 (continued) 86 Block 4 k-I Y Block 3 k-I 1-1-35 (g) Block 3, 4 - k= 1, 12, 18 to show the grid distribution on the leading and trailing edge Figure 5.6 (continued) Block 4 k-18 c . kY Block 3 k 18 Block 4 k-12 Y Block 3 k-12 87 (a) Base grid -- Block 3, 4 -- 40x35 Block 2, 5 - 40x13 x Y z Refined Grid - 79x35 on wing surface (b) Refined grid - Block 3, 4 -- 79x13 Block 2, 5 -- 79x13 Figure 5.7 Base grid and refined grid x Base Grid - 48x35 on wing surface I CHAPTER VI SOLUTION METHOD The multi-block computational method considered is rather simple and straightforward in application. For the wing-fuselage flow problem, the flow field has been properly divided into six contiguous blocks, and the grid network has also been generated. In the solution process, each block is then treated as an independent flow problem that is solved by the same thin-layer Navier-Stokes code. The interface boundary conditions between blocks are updated with a distance weighted interpolation before the computations march for another time step. A computer code has been developed and investigated for the six-block wing-fuselage aerodynamics in which the unsteady 3-D TVD thin-layer NavierStokes code is made to a subroutine subprogram. Since the multi-block computer code has to deal with several blocks of different configurations and their interfaces, the code has to be modular and flexible. Also, the code has to adapt to the architecture of the host computer system. A six-block wing-fuselage Navier-Stokes code (WF6B) was then written for the available computer system, Cray-2 of NAS system at NASA Ames Research Center. The special feature of the supercomputer, Cray-2, is its huge main memory, however, the computing speed is a little slower than Cray Y-MP. In order to take advantage of the large main memory while not sacrifice the computing speed, WF6B uses part of the 88 89 main memory for storing the intermediate data. The other main memory is then used as active area for executing computations. A set of common arrays is then allocated at the active area and used for the computations of each blocks. A block diagram to show this data arrangement is given in Figure 6.1. A detailed description of WF6B can be found in [43]. In order to illustrate the multi-block solution method a flowchart of WF6B is given in Figure 6.2. As shown in the figure, two do-loops run over each blocks and the number of time steps requested. In the inner loop, for each block, BC is used to set up the boundary conditions and then TLNS is to formulate the governing equations and perform the three inversions. A value of L2 residual for each block is computed by RESID every 10 time steps. After the computations run over each block, INTERBC updates the interface boundary conditions, and then the computation continues to the next time step. Finally, the pressure and the shear stress coefficients are computed by CPP and CSP, respectively. In what follows, the techniques used to calculate the metric coefficients and Jacobians of transformation are described first. Then the boundary conditions imposed on the boundary surfaces of the six-block flow problem are discussed while the interface boundary conditions are given in a separate section. Finally, the equations used to compute the L2 residual and pressure and shear stress coefficients are given. 6.1 Jacobian and Metric Coefficient The metric coefficients and Jacobians of transformation can be calculated from Eqs.(2.16-2.17), if the covariant metrics, x,, y,, z,, etc, are known. For a second-order approximation, the covariant metrics are calculated at interior grid 90 Block 1 28x25x60 rl(3,j,k,l) ql(j,k,1,6) Block 2 40x13x60 r2(3,j,k,l) q2(j,k,1,6) Block 3 40x35x60 r3(3,j,k,l) q3(j,k,l,6) Block 4 40x35x60 r4(3,j ,k,l) q4(j,k,1,6) Block 5 40x13x60 r5(3,j,k,l) q5(j,k,1,6) Block 6 26x25x60 r6(3,j,k,l) q6(j ,k,l,6) Note: r's are the Cartesian coordinates q's are the unknown vectors Block data sweep in Active Area 40x35x60 Common Blocks /varl/ /vars/ /btrid/ a(3,j,1,5,5) /vis/ turmu(j,k,l) /var3/ xx(9,1) etc. r-Cartesian coordinates q-unknown vector s-RHS vector a-5x5 matrices turmu-eddy viscosity xx-metric coefficients results sweep out read in grid data and/or restart file write out results Figure 6.1 Data arrangement in WF6B r(3,j,k,1) s(j,k,l,5) q(j,k,1,6) Permanent Disk stores r & q for each block 91 do N = 1,NMAX grid data restart data do BLOCK = 1,6 BC TLNS L E. U INTERBC residual write out r 9 restart data pressure coefficients shear stress coefficients Figure 6.2 The flowchart of WF6B -1 1 H- + - if CPP = yes CPP if CSP = yes CSP read in grid data compute Jacobians & generate initial condition or read in restart data -- I I --- 92 points by a central-difference formula and second-order one-sided differencing at the boundaries, e.g., X = (x,+1,,-xj.1,)/2 for 1< j xe = (-3x,k+4xj+j-xj+2,)/2 for j = 1 (6.1) xf = ( 3xj4-4x._W+xj_.2)/2 for j = jmax The Jacobian, J, of each grid cell is then computed from Eq.(2.17). The contravariant metric coefficients are defined in Eq.(2.16), and can be calculated accordingly. The Jacobians and the metric coefficients evaluated by the above finitedifference approach are rather sensitive to grid characteristics, and can produce a large truncation error if the grid network does not have good properties of smoothness and orthogonality. It is known that they can also be approximated from geometric relations such as cell volume and cell-face area. In the discretized sense, the Jacobian defined in Eq.(2.17) is equivalent to one over the cell volume, 1/Vol, and the contravariant metric coefficients defined in Eq.(2.16) are equivalent to one over the cell-face normal. If one defines a computational cell that is centered at the corresponding grid point, with its cell-faces located half-way to the neighboring points, it is then appropriate to replace the Jacobian, J, by the corresponding cell density, 1/Vol. Furthermore, since ,/J, 4,/J, and c./J defined the x, y, z components of the normals (not unit normals) to the local tangent plane of the constant ( surface, one can replace , , by n/Vol, ny/Vol, and n,,/Vol, where ny, nr and nj are the components of the normals to the local constant j planes. Likewise, one can replace t x, '7, 7. and ., Y, . by 93 nh/Vol, nh/Vol, nh/Vol and nb,/Vol, n1y/Vol, n./Vol, respectively. A general hexahedron, defined by eight arbitrary corner points, can be partitioned into six tetrahedrons, Figure 6.3a. The volume of a tetrahedron denoted by its vertices (ab,c,d) equals one-sixth of the triple product of the three vectors emanating from one of the vertices and ordered according to the righthand directions, Figure 6.3b. For example, V'e(a,b,c,d) = [rd -(rbd x rd)]/6 (6.2) = [(xe-x+(b-Yxa-z -xd)(y-yd)(ze-zd) + (xa-xd)(ye-yd)(zb-zd)-(Xa-Xd)(yb-yd)(Zc-Zd) - (x,-xd(y-Yd)s-zd)- -xXy-yd)(z.zd)]/6 where r, = r - ri, and ri and rj are the position vectors. If the eight vertices of the hexahedron are denoted by numerals 1-8 and are associated with the j,k,l triplets as follows: 1 - j,kl 2 - j+1,k,l 3 - j+1,k+1,1 4 - j,k+1,1 5 - j,kl+1 6 - j+1,kl+1 7 - j+1,k+1,1+1 8 - j,k+1,1+1 the volume of the hexahedron is Vol(j,k,l) = Vet(1,2,3,7) + Vtel(1,2,7,6) + Ve(1,5,6,7) + V'l(1,3,4,7) + Vt(1,4,8,7) + Ve(1,5,7,8) (6.3) The Jacobian of a grid-centered computational cell, which is bounded by the surfaces at the middle of two adjacent grid points, is computed as one over the average volume of all the associated hexahedrons. 94 Vol(j,k,1) = V"(1,2,3,7) + V'(1,2,7,6) + V'(1,5,6,7) + V*0(1,3,4,7) + Vm(1,4,8,7) + Ve(1,5,7,8) (a) Partition of a hexahedron a ~ a--------------b C V"(a,b,c,d) = [r -(rd x rd)]/6 (b) Volume computation of a tetrahedrons Figure 6.3 Volume computation of a hexahedron 5 6 2 3 / 3 |

Full Text |

PAGE 1 NUMERICAL SIMULATION OF TRANSONIC TURBULENT FLOWS OVER A COMPLEX CONFIGURATION SHU-CHENG YANG A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTL^L FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1990 PAGE 2 ACKNOWLEDGEMENTS The author wishes to express his sincere appreciation to his dissertation advisor. Dr. Chen-Chi Hsu, for his guidance and encouragement throughout this endeavor. Special thanks are expressed to Dr. W. Shyy for his invaluable criticism of the draft of this dissertation and his stimulating discussions. Many thanks are also extended to the other supervisory committee members. Dr. B.M. Leadon, Dr. U.H. Kurzweg and Dr. A.K. Varma for their interest and comments on the research project. . . Â» ^ i . : Â• Thanks are also due to Dr. M.H. Chen and Dr. N.H. Shiau for their assistance on their previous work of the TVD scheme and computer codes. Special thanks must go to Mr. J. Ordonez for his tireless proofreading of the manuscript. Also, the author has a debt of gratitude to Dr. N.C. Yao and Dr. C.G. Tu of the Chinese Navy, without whose encouragement and help he would not have had a chance of pursuing the advanced study here at the University of Florida. Finally, the author is indebted to his family for their continuous support. This work was partially supported by a NASA-Ames research grant, NAG2473. Most of the calculations were performed with a Cray-2 of the NAS System at NASA Ames Research Center. ii PAGE 3 TABLE OF CONTENTS Pa ge ACKNOWLEDGEMENTS ii ABSTRACT v CHAPTERS I INTRODUCTION 1 n GOVERNING EQUATIONS 7 2.1 Navier-Stokes Equations and Equation of State 7 2.2 Non-Dimensionalization 12 2.3 Transformed Navier-Stokes Equations 14 2.4 Thin-Layer Approximation 17 2.5 Turbulence Model I9 m NUMERICAL ALGORITHM ?. 25 3.1 Beam and Warming Scheme 26 3.2 TVD Scheme 32 3.3 Flow Solver 43 IV THE FLOW PROBLEM AND SOLUTION STRATEGY 45 4.1 The Flow Problem and Preliminary Work 45 4.2 Solution Strategy 5I 4.3 Six-Block Grid Topology 54 V GRID GENERATION 62 5.1 Surface Grid Generation 63 5.2 Volume Grid Generation 66 5.3 The Six-Block Grid 74 11 PAGE 4 VI SOLUTION METHOD 88 6.1 Jacobian and Metric Coefficient 89 6.2 Boundary Condition 95 6.3 Interface Boundary Condition 105 6.4 Convergence Criterion 108 6.5 Pressure and Shear Stress Coefficient 109 Vn RESULTS AND DISCUSSION Ill 7.1 The Result from Base Grid 112 7.2 The Result from Refined Grid 115 Vm CONCLUDING REMARKS 130 APPENDICES A JACOBIAN MATRICES 132 B EIGEN-FUNCnONS 134 C CLUSTERING FUNCTIONS 136 REFERENCES 139 BIOGRAPfflCAL SKETCH 144 iv PAGE 5 Â•'f 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 TRANSONIC TURBULENT FLOWS OVER A COMPLEX CONFIGURATION By Shu-Cheng Yang December 1990 Chairperson: Chen-Chi Hsu Major Department: Aerospace Engineering, Mechanics and Engineering Science A multi-block computational method is developed for three-dimensional high speed turbulent flows over complex configurations. In this method, the flow domain is divided into several contiguous blocks such that each block is partially bounded by a soUd surface. The chosen grid topology has advantages in computing eddy viscosity as well as in applying a thin-layer Navier-Stokes code. In the solution process, each block is then treated as an independent flow problem governed by the same set of time dependent thin-layer Navier-Stokes equations. The interface boundary conditions between blocks are updated at every time step of the solution algorithm. The application of the method with distance weighted interpolation for updating the interface boundary condition is rather simple and straightforward; however, special measures in grid topology and V PAGE 6 generation are required. The flow solver employed is a total variation diminishing thin-layer Navier-Stokes code implemented with an algebraic eddy viscosity model of Baldwin-Lomax. In this study, the proposed method is investigated for the simulation of a transonic turbulent flow past a wing-fuselage configuration in which the flow field is properly divided into six blocks. A coarser grid was first used for flow simulation, followed by grid refinement smdies. With the use of the refined grid, the calculated wing surface pressure distribution agrees well with measured data. Numerical results obtained also show that the computed shock location and pressure distribution on the wing can be in good agreement with the experiment data if the block grids used are adaptive to important solution characteristics. It is concluded that the proposed method is indeed a very promising approach to be developed further for complex configuration flow simulation. vi PAGE 7 CHAPTER I INTRODUCTION For a high speed flight vehicle, numerous complex three-dimensional fluid flow phenomena appear, including the formation and shedding of vortices, the shock-wave boundary-layer interaction and induced flow separation, and the shock on shock interaction. Since the performance of a flight vehicle is strongly affected by these complex flow phenomena, accurate prediction of them and associated aerodynamic forces and moments are critical for advanced design of a modem flight vehicle. At present an accurate wind-tunnel simulation of the complex flow problem in flight conditions can be very difficult, if not impossible. With the advent of supercomputers and the advancement of solution techniques for nonlinear problems, computational fluid dynamics (CFD) is under extensive development and has been apphed to complement the wind-tunnel experiment and field test in the design process of a flight vehicle. The use of CFD in the aerodynamic design of flight vehicles has / progressed rapidly over the last few years. In the early days, CFD was used to support the vaUdity of a design that developed in the wind-tunnel by trial and error. The state-of-the-art of the CFD method has progressed to the point of being regarded as an important design tool for many flow problems [1]. The reduction of wind-timnel occupancy is only a direct benefit of CFD. Detailed 1 PAGE 8 flow information that is unattainable from the wind-tuimel test often can be obtained by numerical simulations. This information about the underlying flows can lead to a better understanding and ultimately to a better design. In fact, the implementation of CFD has fostered a revolution in the design process of flight vehicles. Current CFD methods have only demonstrated an ability to simulate flows about complex geometries with simple physics or about simple geometries with more complex physics. The panel methods are unique in that they provide a capabihty for solving the flow about completely general configuration. Their principal limitation is that they are restricted to simple physics as modeled by the linear equation. The Euler code provides more flow information than the panel method; however, the important viscous effects are excluded in the approach. For more realistic flow simulations, the application of Navier-Stokes equations is required. However, the use of full Navier-Stokes equations is simply too expensive to be practical for today's computers. An alternative is to use the thinlayer Navier-Stokes equations. In general, the governing equations are approximated by a numerical scheme and then solved in either a structured or an unstructured grid network. The latter is quite flexible in gridding but it requires more work in bookkeeping and computer coding. In the structured grid approach, the governing equations are first approximated by either a finitedifference or a finite-volume scheme, and then solved in a boimdary-fitted curvilinear coordinate system. This approach has advantages in boundary condition treatment and computer coding; however, it may have difficulty in PAGE 9 3 generating good grid system for complex configurations. As evidenced by papers published in journals and presented at conferences, e.g. [2-4], the flow simulations involving simple geometries have progressed considerably; however, the capability for turbulent flow problems involving complex 3-D configurations is still in a development state. Hence, further research and development effort is needed to advance the CFD methods for more complex geometry flow simulations. For flow past a complex configuration, such as a flight vehicle, the generation of an acceptable single structured grid for the entire flow field is rather difficult. It is much easier to divide the flow field into several subregions, termed blocks, and then generate a good grid for each subregion. Hence, it is natural to develop a solution method based on the multi-block grid network. There have been a number of multi-block solution techniques proposed for the Euler simulation of aircraft models and other complex configurations [5-7]. However, the topology of block gridding for viscous turbulent flow simulation is more involved than that for inviscid flows and requires a lot more grid points to resolve the thin viscous layer. In general, a good block grid topology for Euler simulation may not be a proper one for Navier-Stokes sunulation. There are also zonal methods developed to use different equation sets in different regions of the flow field. In general, the flow field is divided into inviscid and viscous zones, and then the Euler equations are appUed in the inviscid zone, while the boundary-layer or Navier-Stokes equations are used in the viscous zone. A zonal approach that was based on the Euler/thin-layer NavierStokes equations was proposed and tested successfully for transonic wing flow [8]. PAGE 10 Another zonal algorithm that applied the Euler and thin-layer Navier-Stokes equations for simulating the flow field of isolated wings was reported in [9]. The flow field was divided into four zones; two inviscid zones with coarse grids and two viscous zones with clustered grids. The computed results were in good agreement with experimental data for the surface pressure, except in the immediate vicinity of the tip and in the shock-induced separated region. The zonal approach also has been used to simulate flow problems for more complex geometry. A transonic viscous flow past an F-16A wing-fuselage configuration has been simulated by a zonal approach [10]. The flow field was divided into as many as 16 zones; in the inner zones adjacent to no-slip surfaces the thin-layer Navier-Stokes equations are solved, while in the outer zones the Euler equations are used. The prediction of wing surface pressures was quite good, except at the leading edge and the aft-shock position. The zonal approach has been shown to be rather efficient. However, the difficulties are to locate properly the interfaces for the inviscid and viscous zones and to patch the solutions that are obtained from different equation sets at these interfaces. Recently, a zonal method with Chimera overset grid scheme was investigated for subsonic turbulent flow about the F-18 fuselage forebody and the combined wing-fuselage [11]. Special coding was required for the overlapped grid region, but the computed results have been shown to be in good agreement with flight-test flow visualization and surface pressure measurements. In this study, a novel multi-block computational method is proposed and explored for the simulation of transonic turbulent flows over complex PAGE 11 configurations. In this method, the flow domain is divided into several contiguous blocks such that each block is partially boimded by a solid surface. For ease in applying thin-layer approximation and the Baldwin-Lomax turbulence model, the sohd smface is mapped onto an entire boundary plane in computational space. In the solution process, each block is then treated as an independent flow problem governed by the same set of time dependent thin-layer Navier-Stokes equations. The interface boimdary conditions between blocks are updated at every time step of the solution algorithm. The solution algorithm with distance weighted interpolation for updating the interface boundary condition and the computer coding are rather simple and straightforward. The multi-block method offers several distinct advantages over the single block computation. In fact, the use of block grids makes the problem of simulating flow fields about complex geometry more tenable. If a geometric feature is changed, only the related block requires to be modified without completely redoing the basic grid topology. Also, the solution procedures are adapted to the modern computer architecmre; in parallel processing, each block can be solved at different processors; or in case * short of main memory, each block can be solved at a time while the other blocks remain on extended storage. A transonic turbulent flow over a wing-fuselage configuration is investigated for the development of the proposed method. The flow field is properly divided into six contiguous blocks. A coarse base grid was first generated and used for the flow simulation, and then a refined grid was generated from the base grid by halving the grid spacing along the wing in the PAGE 12 6 chordwise direction. A flow simulation package that includes grid generation, flow solver, and plotting program, was developed and applied to the six-block flow problem. The method was first investigated for a transonic tiu-bulent flow of Mach 0.8 past the wing-fuselage configuration at zero angle of attack, then simulations for different angles of attack, a =4" and a =-3Â°, were conducted to gain further insight and understanding on the solution algorithm. With the use of the refined grid, the calculated wing surface pressure distribution agrees well with measured data. Nimierical experiments conducted also show that special measures and care are required in generating good grid systems involving excessive distortion between the physical and computational domains; however, the results obtained indicated that the proposed method is a very promising approach for the simulation of complex configuration aerodynamics. PAGE 13 CHAPTER n GOVERNING EQUATIONS Some theoretical background that required for transonic turbulent flow simulations are described in this chapter. The Navier-Stokes equations and ideal gas equation of state are first discussed. Then, they are non-dimensionlized and transformed to a curvilinear coordinate system for ease in numerical appUcations. The thin-layer Navier-Stokes equations valid for high Reynolds number flows are obtained by neglecting the diffusion terms along the du-ection of the soUd surface. The time-averaged thin-layer Navier-stokes equations are used for turbulent flow simulations and the Baldwin-Lomax model [12] is used as the turbulence closure. 2.1 Navier-Stokes Equations and Equation of State 2.1.1 Navier-Stokes Equations The fundamental equations of fluid dynamics are based on the following universal laws of conservation: (a) conservation of mass; (b) conservation of momentum; (c) conservation of energy. The Navier-Stokes equations that govern most of the compressible flow problems can be obtained by applying those universal laws to a fluid flow [13-14]. 7 PAGE 14 8 Continuity equation . The conservation of mass law applied to a fluid passing through an infinitesimal, fixed control volume without sources or sinks of mass yields dp dpyx dps apw ^ , , Â— + + + = 0 (2.1) at ax ay az ^ ' where p is the fluid density and u, v, w are the components of fluid velocity in x, y, z direction, respectively. Momentum equations . The conservation of momentum law appUed to the control volume while neglecting body forces reads iÂ£^4--^(pu^+p)+-^+i^ = zx at ax " dy az ax ay az at ax ay az ax ay az a^w apwu apwv a , , ^ 8tÂ„ dr ar +Â— +Â— + (pw^ + p) = S.+ .L_S.+ "'^a (2.2) at ax ay az ax ay az where the components of the viscous stress tensor ry for a Newtonian fluid are given by the constitutive equation au 3v aw. ^ au 2 au av aw '"xx = A( Â— +Â— +Â— ) + 2mÂ— = Â±^(2^-^.^) ax ay az ax 3 ^ ax ay az ^ " 'ay sz 5x ' ' ay 3 '^^ ay az ax ' _ , 3w 3u av aw 2 aw au av , au av r PAGE 15 av aw , aw au , IS is a Here, p is the thermodynamic pressure and fi is the dynamic viscosity. The fluid assumed to be isotropic and satisfying the Stokes condition, A = -2/3 n, which good approximation for flow problems without relaxation process. For air, the viscosity ^ depends mainly upon the temperature and can be estimated by using an interpolation formula based on Sutherland's theory of viscosity, i.e., , T,L5 TÂ„ + S M = f^.(^r g (2.2b) where denotes the viscosity at the reference temperature TÂ„, and S is a constant which assimies the value 198.6"R [15]. ' ' Energy equation. The conservation of energy law applied to the control volume and assuming no external heat addition results -^[(E+p)u]+ i-[(E+p)v]+ 4-[(E+p)w] = ^ (2.3) 01 ax ay az ax ay az where ! . aT ^, = UrÂ„ + Vr^ + WrÂ„+kÂ— ^ . aT '' . ^ = Ury^ + VTyy + Wry, +kÂ— " , (2.3a) aT = UrÂ„+Vr^ + WrÂ„ +kÂ— , PAGE 16 10 Here, E is the total energy per unit volume, T is the temperature and k is the thermal conductivity. The heat loss by conduction qj is related to the temperature gradient by Fourier's law of conduction aT = -kÂ— (2.4) The equations described above, Eqs.(2. 1-2.3), constitute a complete set of the time dependent Navier-Stokes equations that is valid for a class of compressible flow problems without heat addition and chemical reaction (or relaxation). The Navier-Stokes equations includes one continuity equation, three momentum equations and one energy equation, but there are seven unknowns, namely, p, u, v, w, p, E, and T. The relationships between the thermodynamic variables (p, p, T, E) are established by the equation of state of the working fluid. 2.1.2 Equation of State In this study, air is the assumed working fluid and the perfect gas equation of state is employed. For a perfect gas, the thermal and caloric equations of state are p = pRT, e = QT ( Q = constant ) (2.5) respectively. Here, R is the universal gas constant, e is the internal energy. Some other useful relations are h = CpT, 7 = Cp/Q a^ = -yRT (2.6) where h is the enthalpy, Cp and Q are the specific heat at constant pressure and volume, respectively, 7 is the ratio of specific heats, and a is the local speed of sound. PAGE 17 If only internal energy e and kinetic energy are considered significant, the total energy E can be written as E = p[e + (u^+v^+w^)/2] Consequently, one obtains two specific equations relating p, T to the dependent variables p, u, v, w, and E, i.e., 1 c 1 (2.7) If expresses heat conduction qj in terms of internal energy through the use of caloric equation of state, one can rewrite Eq.(2.4) as follows. , aT k ae yu ae Here, the Prandtl number Pr is defined as Pr = C^ti/k, for air Pr = 0.72, approximately, and the internal energy is given by e = ^-yCu'+v^+w^) (2.8a) and the /9's in Eq.(2.3a) can also be rewritten as iii de = UrÂ„ + VrÂ„ + WrÂ„ + " Pr ax ae = Ur^ + Vr^+Wr^ + Â— Â— (2.9) = UrÂ„ + VrÂ„ + Wr + Pr ay 7/i ae " Â° Pr az PAGE 18 2.1,3 Vector Form 12 The Navier-Stokes equations, Eqs.(2. 1-2.3), form a system of five equations and can be cast in vector form as aq aE aF aG aEL aF, aG, at ax ay az ax ay az (2.10) where q is an unknown vector, E, F, G are the flux vectors and E^ F^ G^ are the viscous flux vectors. The explicit expressions for these vectors are q = p \ 1 pU pu +p ' pV E = < pUV , F = pW puw [e J L(E+p)uJ v pV pVU pv^+p pVW G = pW pWU pWV 3W^ + pW^ + P Ue+p)wJ 0 xz I9x J F = '0 " ro (2.10a) ' T zy [fir) in which ry are given by Eq.(2.2a), ^'s are defined in Eq.(2.3a), and p, T are related to q by Eq.(2.7). 2.2 Non-dimensionalization To reduce the number of parameters that characterize the same class of problems, a dimensional analysis is employed. The non-dimensional quantities are obtained by choosing characteristic scales in the following manner: PAGE 19 where the non-dimensional variables are denoted by an asterisk *, the freestream conditions by Â«, aÂ„ is the freestream speed of sound, and L is a characteristic length scale. " ' By substituting the above relations into the governing equations, Eq.(2.10), the non-dimensional equations in vector form are aq* , aE* aF* aG* Â„ dE* dF* aG ' -T+ Â— s+ Â— s+ Â— r = Re X Â— ^ + Â— ^ + Â— /) (2.12) at ax ay az ^ ax ay az ^ . ; ^ ^ where the Reynolds number, Re, is defined as v. Re = and the relationships from the perfect gas equation of state, Eq,(2.7) become p* = (^-l)[E*Ap\u*^+v-2+w-2)]= (-y-l)pV C V Â• E* 1 ' " 1^ = 7(7-l)[Â— -^(u'^+v'^+w'^)] ''^^ ^ Â• (2.13) The unknown vector q* and flux vectors E*, F*, G*, E/, F/ and G* have the same form as in Eq.(2.10a), except that all the quantities are now non-dimensional. In summary, the dimensionless equations govern a class of flow problem. Only two parameters, namely, the Reynolds number, Re, and the Prandtl number, Pr, appear in the equations, the Mach number is dropped out due to the use of a, as characteristic velocity scale. In the following, the asterisks * are removed from the non-dimensional equations except where noted. PAGE 20 14 2.3 Transformed Navier-Stokes Equations In the finite-difference approach, to enhance numerical accvu-acy and efficiency in applying boundary conditions, the governing equations are transformed from the Cartesian coordinates to a boundary-fitted curvilinear coordinate system. Consequently, a general computational code can be developed such that the computation is performed in an uniformly-discretized computational domain [16]. If one introduces a generalized coordinate system e = ?( x,y,z,t ) n = r,( x,y,z,t ) (2.14) r = r( x,y,z,t ) r = t and applies the chain rule of partial differentiation, the non-dimensional governing equation, Eq.(2.12), becomes q. + QeC. + q.-^. + q^ r . + (C^^ + v^fi, + ^^,) + (CyF^ + VyF, + TyF,) + (e,G^ + r,fi^ + r,G,) = Re-^[(^ A, + r,^^ + rxE^) + + VyF^ + CyF^) + (C,G^ + vfi^ + r,G^)] (2. 15) The metric coefficients appearing in these equations can be determined by the matrix relation as follow i^ ' ^? ^f) Vy Â»?, y? y. Vc yr k '<} . T.f ZÂ„ Zf. . o' 0" \ . or PAGE 21 15 f/x 'y X ''y ''z fx f y f z . Xj. y? y. yf Zf -1 = J (yfZf-yfZe) /fZ^-x^z -(x^y^-x,y^) (XfZ^-x,Zj) x^y^-x^y^ J yez.-y.Zf '^x Cy ^0 ' Â»Â»Â» '?x '?y '?Z "yr ' I Tx Ty fz J Xj X^ Xj. y? y. yf z. where J is Jacobian of the transformation and J'\ representing the volume of the grid cell in physical space, is given by _ a(x,y,z) a(e,i7.f) = Xe(y.ZryfZ,)-x,(yeZf-yfZ^)+x^(y^z^-y^Zj) . . (2.17) In order to put the transformed equations, Eq.(2,15), into conservation-law form, the equation is first divided by the Jacobian and then rearranged to give (qJ/J + (l.qf + ?xEe + eyF, + e.G,)/J + (Â„,q^ + ^ JE, + + ^^g,)/j + (r.qf + rxE, + fyF, + f ,G,)/J = Re-^[(e + e/,, + ?,G,,) and then rewritten into the following conservation law form (q/J)r + (a.q+e,E+eyF+C,G)/J)^ + (M+Â»7^+r,yF+r,,G)/J)Â„ + ((r.q+fxE+fyF+r,G)/J), = Re-^[((eA+^yF,+ ?,GJ/J), + . + ((.7A+'7yF,+ r,,G,)/J), + ((rJE,+ ryF,+ f,G,)/J),] , (2.18) in which the follov^ing metric identities [10] have been apphed. PAGE 22 16 (2.19) Â— (7^)+Â— (f) = 0 a? ^ J ^ at; M ar M ^ The transformed equation, Eq,(2.18), has the same form as the original equation, Eq.(2. 12), and can be written as v V ' .v ar ae aÂ»7 ar a$ a?; ac (2.20) where q = q/J, E = (^,q+|^+C/+|,G)/J, F = (,7.q+Â»7^ + r7yF+r,,G)/J, A G = (r,q+fJE+fyF+f,G)/J, (2.20a) A = (^A+?,Ft+?zG.)/J, A = (Â„A+'/yFT+'7zG,)/J, = (rA+ryF,+f,G,)/j The explicit expressions for the transformed unknown vector, flux vectors and viscous flux vectors are T-i q = J p pV pVi E F = J-^ E = J -1 /juU+^j) pvU+^p pwU+^j) L(E+p)U-^j5j puV+r/J) pvV+Â»7yP pwV+r/J) (E+p)V-Â„j)J 0 Cx'-xy+Vyy + ^z''xy ^x'-xt+Vp'^^^^^ ^x^x+?A+^A G = J-^ puW+fj) pvW+Typ pwW+fj) [(E+p)w-rj)J F, = J-^ 0 '?x'"xx+Vyx+''z'"zx '7x^xy+Vyy + ''z'-zy (2.20b) PAGE 23 0 where U, V and W are contravariant velocity components given by V = f7,+T7,u+r;yV+Â»7,w (2.20c) w = r,+f,u+fyV+r^w ' It is noted that the velocity gradients in Ty, Eq. (2.2a), and the temperature gradient in /9's, Eq.(2.3a), have to be transformed into computational space also, e.g., Â• -k ax ay az ' ax ^r^ 3u au au av^ av av ^ a? aÂ»; ar . aw ^ _ aw aw ^ au au au + (Â— Â— Â»7z+ Â— rz)]+24(Â— -Ffx)] ae aÂ»7 ar a^ at? ar 2.4 Thin-Layer Approximation The numerical simulation of complete Navier-Stokes equations, Eq.(2.20), is in general rather time consuming and needs a huge amount of computer storage. For high Reynolds number flows, however, the viscous effects are often confined to a thin layer near no-slip boundaries, called the boundary-layer, outside of which the flow can be considered to be mostly inviscid. With this in mind, it is reasonable to assume that the diffusion processes along the body PAGE 24 18 surface are negligibly small in comparison with those in the direction normal to the body surface. 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, while all other terms in the momentum equations are retained. One of the advantages of retaining the terms which are normally neglected in boundarylayer theory is that separated and reverse flow regions can also be computed. Also, flows which contain a large normal pressure gradient can be readily computed. Another justification of the use of a thin-layer approximation can be given from a computational viewpoint. In practice, a substantial fraction of the available computer storage and time is expended in resolving the cross-stream gradients in the boundary layer since a highly stretched grid is required, and the resolution along the body is treated as what is needed in inviscid flow. As a result, even though the full derivatives are retained in the equations, the gradients along the body are not resolved in an adequate manner. Hence, one can drop those terms that are not being adequately resolved, provided that they are reasonably small. Since the thin-layer model requires a boundary-layer type of coordinate system, it is assumed that the sohd body surface is mapped onto an entire CÂ»7plane in the transformed space, then the thin-layer approximation requires that the e and f; derivatives in the viscous terms be dropped. Upon simplifying the complete Navier-Stokes equations, Eq.(2.20), the system of governing equations becomes i PAGE 25 aq , aE , aF , aG _ aS dr where ae dri = Re-\Â— ) ac (2.21) S = J-^ 0 * miUj. + m25-, miVf + mzfy I mim3 + m2(r^u + ryV+r,w) J (2.21a) ^2 = M(r,Uf + fyVf + r,w^)/3 (2.21b) mj = (u2+v2+w2)^/2 + PT\y-ir\r), Note that the notation has been removed from the above equations for simplicity and the heat conduction term are replaced by jfi ae _ /i aT '^""Pr ar ~ "(7-l)Pr "aT ^^'^^^ The thin-layer Navier-Stokes equations, Eq.(2.21), in conjunction with the equation of state, Eq.(2.13), constitute the required governing equations for this study. They are valid for both laminar and turbulent flow problems. However, due to the widely disparate scales of turbulence, it is impractical to apply the equations directly for turbulent flow simulation. An alternative is to use the time-averaged equations coupled with a turbulence model. 2.5 Turbulence Model The time-averaged equation is obtained by decomposing the dependent variables into time mean and fluctuating components and then time-averaging the PAGE 26 20 entire equation. In the averaging process, new unknowns will be introduced and must be treated by some approximating techniques so as to close the problem. In this study, the mass-weighted time-averaged equations are used for turbulence flow simulation with the eddy viscosity that appears in the averaged equations provided by the Baldwin-Lomax two-layer algebraic model. In this approch, the mass-weighted time-averaged thin-layer Navier-Stokes equations will have the same form as the original equations, Eq.(2.21), while all the dependent variables are time-averaged. The effective coefficient of viscosity n is comprised of two parts, the laminar viscosity /i, and the tiu-bulent eddy viscosity /i,. The laminar viscosity A*, is obtained from Sutherland's theory of viscosity, Eq.(2.2b), and the eddy viscosity ^l^ is provided by a turbulence model [16]. Similarly, the effective thermal conductivity k is also comprised of two parts, k, and k,. The heat conduction is expressed in dimensional quantities as The heat conduction term in Eq.(2.22) is then replaced by its dimensionless counterpart Here, the turbulent Prandtl number Pr, is defined as Pr, = ^x^c^/k^, and is about 0.9 for air [10]. In this approach, the eddy viscosity /i, is the only quantity that has to be modeled. In what follows, the theoretical background about the algebraic model will be discussed first, and then the Baldwin Lomax model is described. PAGE 27 21 2.5.1 Two layer algebraic model In analogy with the coefficient of viscosity in Stoke's law for laminar flow, Boussinesq [17] introduced a mixing coefficient, /iÂ„ or referred to as eddyviscosity, for the Reynolds-stress that appears in the time-averaged equation, au ^ = /^.-TÂ— (2.24) ay where u is the streamwise velocity, y is the spatial coordinate normal to the solid siuface, and 3u/ay is assumed to be the only significant component of strain rate. With Prandtl's mixing length hypothesis, the coefficient n^, is then related to mean velocity field by 'vdu dy . . . :r (2.25) where 1 is the mixing length. ' The algebraic model is then to establish a relation between the mixing length and the characteristic length of the respective flow. For simple wall shear flows, the mixing length in the law of the wall region is proportional to the normal distance from the wall, i.e., 1 = ky, where k = 0.41 is a universal constant determined experimentally. In the viscous sublayer, however, the mixing length theory is not valid. A corrected form of the mixing length valid down to the wall was deduced by van Driest. He introduced a damping factor and wrote 1 = ky[l-exp(-yVA^)] . , (2.26) PAGE 28 22 where y"*" = ypujn is the law of the wall coordinate, A"^ denotes an empirical constant which equals 26 for flat-plate flow, and u^, the friction velocity, is defined as = {rjpf', where rÂ„ is the shear stress on the wall. In the law of the wake region, the Clauser model is commonly used in conjunction with a factor to account for the intermittent effect. The eddy viscosity is then written as where the value is obtamed experimentally and assumed to be 0.0168, U is the streamvdse velocity on the edge of the boundary layer, 5* is the boundary layer thickness, 7 is the intermittency factor deduced from a curve fit of measured data and given by Reynolds and Cebeci [18] as where s is the boimdary layer thickness. Based on the above arguments, Cebeci and Smith [19] proposed a well known two-layer algebraic model as follows. In the iimer-layer: (2.27) 7 = [l + 5.5(y/5)'] (2.27a) (2.28) 1 = ky[l-exp(-yVA*)] (228a) In the outer-layer: A*, = k,U5*7 (2.28b) 7 = [l+5.5(yA)Â«] (2.28c) PAGE 29 This model has been shown to be in excellent agreement with experimental ' ; , observations for simple shear flows. However, the model requires the boundary layer thickness and the edge velocity which constitute a practical disadvantage in the implementation of the model. In an attempt to overcome this difficulty, Baldwin and Lomax [12] modified the outer layer formulation to eliminate the necessity of finding the boundary layer thickness. In addition, they replaced the velocity gradient by the vorticity, which is invariant under coordinate transformation, so that the model can be easily applied to 3-D problems. The Baldwin-Lomax model is described next, 2.5.2 Baldwin-Lomax Turbulence model The Baldwin-Lomax model is similar to the Cebeci-Smith two-layer algebraic turbulence model. The eddy viscosity in a profile is calculated for each inner and outer layer, and then matching both layer as follows [12]: where y, is the smallest value of y at which values the inner and outer formulas are equal. The inner layer follows the Prandtl-Van Driest formulation where y is the normal distance to the surface and |w| is the magnimde of the vorticity (2.29) (/^.)inÂ»er = pl'kl with 1 = ky[l-exp(-y VA*)] (2.30) (2.31) PAGE 30 24 = y(py,^J^/f^ is the law of the wall coordinate, A"^ = 26, and /x^ are the local density, shear stress, and laminar viscosity evaluated at wall, respectively. The eddy viscosity for the outer layer is given by (M,)outer = pk,QpIVakeFkleb(y) ' ' ' Â• ' " (2.32) where = 0.0168 is the Clauser constant and = 1.6. * . H ' , The parameter F^^^^ is given by Fwak. = mini J^rf /p where = 0.25 (232a) and Udif = (u'+v^+w^)^^ (u^+v^+w^)^^ . (2.32b) The quantities y^ and are determined from the function F(y) = yk|[l-exp(-yVA*)] (2.32c) where is the maximum value of F(y) and y^ is the value of y at which F^ occurs. The Klebanoff intermittency factor Fu,,,(y) is given by Fkid.(y) = [l + 5.5(Cu^/y^']-Â» and C^.^ = 0.3 (2.32d) The Baldwin-Lomax model has been investigated and apphed to a variety of 2-D and 3-D flow calculations with satisfactory results [10]. The model is capable of handling the compressible flow problems with modest flow separation [20]. The implementation of algebraic models requires only a small amount of computer time and storage; this is particularly important in 3-D computations. PAGE 31 CHAPTER in NUMERICAL ALGORITHM In numerical simulations, the transformed governing equations described in the previous chapter are approximated by finite-difference equations and then solved in a equally spaced computational space. The basic computer code used in this study is a TVD thin-layer Navier-Stokes code which is based on an original version of thin-layer Navier-Stokes code, ARC3D of NASA Ames Research Center. The original ARC3D was based on the Beam and Warming scheme [21], and has been improved for its efficiency, accuracy and convergence [22]; for instance, there are options for using different numerical dissipations for better stability criteria. Alternatively, to improve the robustness of the code, a symmetric Total Variation Diminishmg (TVD) scheme was implemented into the original ARC3D. This modified ARC3D has shown to be more efficient and robust than the original one [23-26], and hence is adopted in this study. _ The BeamWarming scheme and the TVD scheme were constructed by totally different philosophies; however, in applications they can be viewed as the same scheme with different numerical dissipation terms [27]. Both schemes will be described and discussed in the following. ' 25 PAGE 32 26 3.1 Beam and Warming Scheme The BeamWanning scheme is an non-iterative, implicit approximate factorization finite-difference scheme. The implicit operator is factorized (spatially split) to obtain an ADI-like ( Alternating Direction ImpUcit ) algorithm such that the three-dimensional inversion process is reduced to three onedimensional processes. The basic algorithm is second-order accurate in both time and space [28]. The implementation of this scheme to the system of governing equations, Eq.(2.21), is described as follows. 3.1.1 Temporal-Differencing To approximate the time derivative of Eq.(2.21), the well-known secondorder Trapezoidal Rule or Crank-Nicolson Formula is employed to give ^ = + 3-7^) + 0[(Ar)2] where Aq-=q-*V (3.1) If one rewrites Eq.(2.21) as and then substitutes it into Eq.(3.1), the following equation results 2 " 3{ 3, ac ^ 3f In the above equation, the flux vectors E, F, G, and viscous flux vector S at the advanced n+ 1 time step are in general nonlinear functions of the unknown PAGE 33 ; if 27 4 vector q. They are linearlized while maintaining second-order time accuracy as follows. The flux vectors are linearlized by using a Taylor series expansion about q" as follows E"+^ = E"+A"Aq" + 0[(Aq)^] where pn+i ^ pn^gn^qn ^ o[(Aq)^] wherc GÂ°** = G"+C"Aq" + 0[(Aq)^] where The above approximation is second-order in time since Aq" = q"-'Â»-q" = 0(Ar) and 0[(Aqf] = 0[(Arf] The viscous flux vector S is linearized by using the procedure suggested by Steger [29]. The coefficient of viscosity ^ and thermal conductivity k are assumed to be locally independent of Jq ( recall that Jq is the unknown vector in physical domain ), so that elements of S can be expressed in the general form Sj = J'^aj -Â— i where a, is independent of Jq. Each element is then linearized in the following manner by Taylor series expansion S,-^ = S,Vj-l{a,"-^[E(-^r]}jAq/ + 0[(Ar)^] Written in vector form, the viscous flux vector becomes gn+i ^ s"+J-^M"JAq" + 0[(Ar)2] (3,4) A = B = C = 3E aq dF aq aG aq (3.3) PAGE 34 28 The matrices A, B, C and M are so called Jacobian matrices and the elements of these matrices are given in Appendix A. After substituting all the linearized terms, Eqs.(3.33.4) into Eq.(3.2), one has 2 ^ 3^ dr, af ar ^ ^ , ^ . r ... ^ = .^,(iE^_aF^_aG_j^^.,^ , ^ v a^ ar; ar af ^ ^ ^ ^ where I is the identity matrix, and Aq" is the delta form unknown to be solved for. Note also that a first-order time-accurate Euler implicit scheme can be obtamed by simply replacing Ar/2 on the left-hand side of Eq.(3.5) by Ar [28]. 3.1.2 Approximate Factorization Solving Eq.(3.5) directly is very inefficient since it involves the inversion of a very large system resulted from three dimensions. Beam and Warming appUed an approximate factorization procedure to reduce the inversion to three onedimensional inversions. They rewrote Eq.(3.5) in the factored form Since the error introduced by this factorization procedure has the leading thirdorder term Ar^._aA_ _aB_^ aB_ _aC aC aA ag" 4 ac dr, dr, ar ar ac the second-order accuracy of the difference approximation is not affected. PAGE 35 3.1.3 Spatial differencing The spatial derivatives appearing on the left-hand side of Eq.(3.6) are approximated by second-order finite-difference formula, e.g. hence the resulting scheme possesses a block-tridiagonal structure that can be solved efficiently. On the right-hand side, the same second-order differencing is employed for the viscous term, while a fourth-order finite differencing approximation, e.g. is used for the convection terms to keep the convection truncation error from exceeding the magnitude of that of the viscous term, ^ ; After applying spatial differencing into Eq.(3.6), the spatial derivatives are replaced by difference operators, and the finite-difference form of Eq.(2.21) is then . .. " , [1+ y-6,A]"[I+ y^5,B]"[I+ ^(5^C-Re-^S,(J-^MJ))]"Aq= -ArC^^E+sZ+SfG-Re-^S)" (3.7) This is the basic Beam-Warming algorithm. The implicit operator has been decoupled to produce an ADI-like scheme, and due to the second-order central-difference operators employed, the algorithm produces block tridiagonal systems for each spatial direction. PAGE 36 30 3.1.4 Numerical Dissipation The BeamWanning scheme requires the use of numerical dissipations to damp out the spurious oscillations that occurs when dealing with discontinuities. Several different numerical dissipations have been investigated and discussed in [30]. Those employed in the original ARC3D code were proposed by Beam and Warming [28], and Steger [29] and are described as follows: The second-order implicit dissipation operators = -Â«iJ"\a^J (3.8) are added to the imphcit operators on the left-hand side of Eq.(3.7) in the |, rj and f directions, respectively, and on the right-hand side the following explicit fourth-order dissipation term is added: De = -^EJ-'[(VA)'+(V/+(VfA,)']JAq(3.9) in which symbols a and v are the forward and backward difference operators, e.g. ^ = + = -J^i^rki) + 0[A|] and Â€j, eg are constant coefficients used to control the amount of numerical dissipations. It was suggested to set Â€E=At and ei=2Â€E so that as the time step increases, the numerical dissipations added relative to the spatial derivatives of convection and diffusion remains constant. PAGE 37 After inserting numerical dissipations into Eq.(3.7), the finite-difference representation of Eq.(2.21) becomes [1+ |l-5,A+D,r[I+ ^6^B + DJ"[I+ y^(5,C-Re-^5,(J-^MJ)) + D,]"Aq" = -Ar(SjE+6^F+5j.G Re"^5j.S)"+DE (3.10) 3.1.5 Solution procedure For simplicity, the following notation is introduced h = [I+y-5,A+D,r LÂ„ = [1+ ^Â«Â„B + D 1" Lf = [1+ y^(5fC-Re-^5f(J-^MJ)) + Df]"Aq" R" = -Ar(5^E+5,,F+5^G-Re-^5^S)"+DE where the L's are operators, then the Beam-Warming scheme, Eq.(3.10), can be written as L,L,L,Aq= R(3,11) If one further defines intermediate solutions as Aq* = Lj-Aq" Aq = L,Aq the unknown vector q at the n+1 tune level can be computed as follows. LfAq= RÂ» L. Â• /q = Aq Vq" = Aq* (3 12) PAGE 38 q"+i = q-+Aq" The Beam-Warming scheme has been shown to be accurate and efficient for simple flow problems [10], however, for more complex flow cases it can be very sensitive to grid resolutions [25]. . : 3.2 TVD Scheme It is known that finite-difference schemes may produce oscillations when appUed across a discontinuity. For instance, the second-order Lax-Wendroff scheme u;*i = -K uj^i-u? )-V2uil-u){ uj"^i-2uj"+u/!i ), u = aAt/Ax (3.13) for a linear hyperbolic equation au + a =0, a > 0. at ax can produce spurious oscillations near the discontinuity. The Total Variation Diminishing (TVD) scheme is one of the methods designed to eliminate those undesirable oscillations. The notion of TVD was first introduced by Harten [31]. He derived a set of sufficient conditions from which the constructed schemes will not generate oscillations across discontinuities. Since the theoretical background of TVD schemes is rather complex and involved, only relevant basic definitions are given here. The construction of a symmetric TVD scheme and its application to thin-layer Navier-Stokes equations is discussed in the subsequent sections. PAGE 39 A numerical scheme, either explicit or implicit, can be written in the following operator form [32] L uÂ°** = R u" (3.14) where u""^* is the mesh function at the advanced time step to be solved, u" is the known mesh function at current time step, and L and R are some finitedifference operators designed for the particular scheme, e.g., for the second-order Lax-Wendroff scheme given in Eq.(3.13) L = 1 R = 1 -.^Vm-MMWh ' V:l Here and in the following, Aj+j^ stands for the central difference operator, i.e., Aj^.j4U = Uj+i Uj. The total variation of a mesh function uÂ° is defined to be =iJL I^Vi-u"! =jL ''^ > ' ' (3.15) Here, the notation Â« is used to emphasize that the summation is taken over the entire domain. If a numerical scheme applies to a mesh function u, the total variation of the mesh function is non-increasing as a function of time, i.e. TV(u"*M PAGE 40 34 It has been proved that a TVD scheme is monotonicity preserving [32], that is, for any monotone mesh function u that increasing (or decreasing) monotonously with time, w = Qu is also monotone, where Q is a finite-difference operator. Since a monotonicity preserving scheme will not generate spurious oscillations, and neither will the TVD scheme. In what follows, the construction of the symmetric TVD scheme for a onedimension hyperbolic equation will be described first, then the scheme is applied to the Euler equations, and finally it is extended to the thin-layer Navier-Stokes equations. 3.2.1 Symmetric TVD Scheme The symmetric TVD scheme was developed by Yee [27] and it is a generalization of Roe [33], Davis [34] and Sweb/s [35] TVD Lax-Wendroff scheme. To illustrate the construction of this scheme, a one-dhnension linear hyperbolic equation of the form 8u du + a Â— at ax + aÂ— = 0 (3.17) is considered. Sweby has shown that the second-order Lax-Wendroff scheme for this equation can be constructed by adding an anti-diffusive flux to the first-order upwind scheme [35], i.e., Uj-** = i.A..j,u" y2i.(l-..)Aj.j4Aj^^,UÂ» , if a > 0. uj*Â» = ^Aj^^u" ViÂ«.(l + ..)A.^^Aj.j,u, if a < 0. (3.18) PAGE 41 where v = aAt/Ax is the Courant number. Then he introduced a limiting function p and rewrote the above scheme as i^*i=uJ-.{l + y2(l-.)[p(r/)/r/-p(rV,)]}Aj.,^^ , if a > 0 ^^^^^ 1^*1= uJ-u{l-y2(l + ..)[p(rj;i)-p(r-p/rj]}Aj,^u, if a < 0 where p is chosen to be a function of the ratio of consecutive cell gradients rj; pj = P(rj). = ^-mu/Aj+mU, rj' = Aj+^^u/Aj.^^u, with the restrictions 0 I 0 if a < ( KYr-) ^^-^ (3.20a) \'^(1 + u)\p(t: )-1] if a < 0 To avoid determining the upstream direction, Davis then modified the dissipation term as follows, = M(i-IH)[i-P(V)]/2 K-(ij-) = IH(l-|H)[l-p(ri-)]/2 . (3.20b) Shortly after that. Roe [33] further generalized Davis' scheme, Eq.(3.20), and suggested a new limiting function Q. Roe's scheme takes the form |