**Suggested Citation:**"Developing an Accurate and Efficient Method for Viscous Compressible Flow Simulations - An Example of CFD in Aeronautics." National Research Council. 1990.

*The Proceedings: Fifth International Conference on Numerical Ship Hydrodynamics*. Washington, DC: The National Academies Press. doi: 10.17226/1604.

**Suggested Citation:**"Developing an Accurate and Efficient Method for Viscous Compressible Flow Simulations - An Example of CFD in Aeronautics." National Research Council. 1990.

*The Proceedings: Fifth International Conference on Numerical Ship Hydrodynamics*. Washington, DC: The National Academies Press. doi: 10.17226/1604.

**Suggested Citation:**"Developing an Accurate and Efficient Method for Viscous Compressible Flow Simulations - An Example of CFD in Aeronautics." National Research Council. 1990.

*The Proceedings: Fifth International Conference on Numerical Ship Hydrodynamics*. Washington, DC: The National Academies Press. doi: 10.17226/1604.

**Suggested Citation:**"Developing an Accurate and Efficient Method for Viscous Compressible Flow Simulations - An Example of CFD in Aeronautics." National Research Council. 1990.

*The Proceedings: Fifth International Conference on Numerical Ship Hydrodynamics*. Washington, DC: The National Academies Press. doi: 10.17226/1604.

**Suggested Citation:**"Developing an Accurate and Efficient Method for Viscous Compressible Flow Simulations - An Example of CFD in Aeronautics." National Research Council. 1990.

*The Proceedings: Fifth International Conference on Numerical Ship Hydrodynamics*. Washington, DC: The National Academies Press. doi: 10.17226/1604.

**Suggested Citation:**"Developing an Accurate and Efficient Method for Viscous Compressible Flow Simulations - An Example of CFD in Aeronautics." National Research Council. 1990.

*The Proceedings: Fifth International Conference on Numerical Ship Hydrodynamics*. Washington, DC: The National Academies Press. doi: 10.17226/1604.

**Suggested Citation:**"Developing an Accurate and Efficient Method for Viscous Compressible Flow Simulations - An Example of CFD in Aeronautics." National Research Council. 1990.

*The Proceedings: Fifth International Conference on Numerical Ship Hydrodynamics*. Washington, DC: The National Academies Press. doi: 10.17226/1604.

**Suggested Citation:**"Developing an Accurate and Efficient Method for Viscous Compressible Flow Simulations - An Example of CFD in Aeronautics." National Research Council. 1990.

*The Proceedings: Fifth International Conference on Numerical Ship Hydrodynamics*. Washington, DC: The National Academies Press. doi: 10.17226/1604.

**Suggested Citation:**"Developing an Accurate and Efficient Method for Viscous Compressible Flow Simulations - An Example of CFD in Aeronautics." National Research Council. 1990.

*The Proceedings: Fifth International Conference on Numerical Ship Hydrodynamics*. Washington, DC: The National Academies Press. doi: 10.17226/1604.

**Suggested Citation:**"Developing an Accurate and Efficient Method for Viscous Compressible Flow Simulations - An Example of CFD in Aeronautics." National Research Council. 1990.

*The Proceedings: Fifth International Conference on Numerical Ship Hydrodynamics*. Washington, DC: The National Academies Press. doi: 10.17226/1604.

**Suggested Citation:**"Developing an Accurate and Efficient Method for Viscous Compressible Flow Simulations - An Example of CFD in Aeronautics." National Research Council. 1990.

*The Proceedings: Fifth International Conference on Numerical Ship Hydrodynamics*. Washington, DC: The National Academies Press. doi: 10.17226/1604.

**Suggested Citation:**"Developing an Accurate and Efficient Method for Viscous Compressible Flow Simulations - An Example of CFD in Aeronautics." National Research Council. 1990.

*The Proceedings: Fifth International Conference on Numerical Ship Hydrodynamics*. Washington, DC: The National Academies Press. doi: 10.17226/1604.

**Suggested Citation:**"Developing an Accurate and Efficient Method for Viscous Compressible Flow Simulations - An Example of CFD in Aeronautics." National Research Council. 1990.

*The Proceedings: Fifth International Conference on Numerical Ship Hydrodynamics*. Washington, DC: The National Academies Press. doi: 10.17226/1604.

**Suggested Citation:**"Developing an Accurate and Efficient Method for Viscous Compressible Flow Simulations - An Example of CFD in Aeronautics." National Research Council. 1990.

*The Proceedings: Fifth International Conference on Numerical Ship Hydrodynamics*. Washington, DC: The National Academies Press. doi: 10.17226/1604.

**Suggested Citation:**"Developing an Accurate and Efficient Method for Viscous Compressible Flow Simulations - An Example of CFD in Aeronautics." National Research Council. 1990.

*The Proceedings: Fifth International Conference on Numerical Ship Hydrodynamics*. Washington, DC: The National Academies Press. doi: 10.17226/1604.

**Suggested Citation:**"Developing an Accurate and Efficient Method for Viscous Compressible Flow Simulations - An Example of CFD in Aeronautics." National Research Council. 1990.

*The Proceedings: Fifth International Conference on Numerical Ship Hydrodynamics*. Washington, DC: The National Academies Press. doi: 10.17226/1604.

**Suggested Citation:**"Developing an Accurate and Efficient Method for Viscous Compressible Flow Simulations - An Example of CFD in Aeronautics." National Research Council. 1990.

*The Proceedings: Fifth International Conference on Numerical Ship Hydrodynamics*. Washington, DC: The National Academies Press. doi: 10.17226/1604.

**Suggested Citation:**"Developing an Accurate and Efficient Method for Viscous Compressible Flow Simulations - An Example of CFD in Aeronautics." National Research Council. 1990.

*The Proceedings: Fifth International Conference on Numerical Ship Hydrodynamics*. Washington, DC: The National Academies Press. doi: 10.17226/1604.

Below is the uncorrected machine-read text of this chapter, intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text of each book. Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.

Developing an Accurate and Efficient Method for Compressible Flow Simulations - An Example of CFD in Aeronautics - K. Fujii The Institute of Space and Astronautical Science Sagamihara, Japan Abstract Capability of the current CFD technology in the aeronautical society is discussed with one of the rep- resentative Navier-Stokes codes in Japan, as an exam- ple. The code is named LANS3D and was developed for the numerical simulation of high-Reynolds number compressible flows. The algorithm used in this code is briefly described First and then the applications for aircraft simulation, cortical flow simulation and spacer plane simulation are presented. These example show how the code has been improved to satisfy two impor- ta~t requirements of the computational fluid dynamics (CFD) codes; efficiency and accuracy. Importance of developing a supporting system on graphic worksta- tions is kin ally discussed. 1. Introduction Computational Fluid Dynamics (CFD) is begin- ning to play an important role in the aeronautical industry all over the world. People now realize that CFD can be a new and effective design tool with the aid of supercomputer. At the same time, CFD is be- coming an important scientific tool for the fluid dy- namics research. As high-speed supercomputers with large memory become popular, the research now is focused on the Navier-Stokes simulations for under- standing fluid physics as well as for the future use as a design tool in engineering. In the application of the CFD to the aeronautical problems, there are two important remarks to note. First, representative Reynolds number is large com- pared to many of the other CFD applications. Even though viscous equations are solved, contribution of viscous terms is small in most of the flow field. Thus, computational algorithms that have been developed for the Euler equations are mainly used. Physically, viscous effect is confined to the thin layers near the body surface and to evaluate viscous effect properly with the Navier-Stokes simulations, grid spacing near the body surface should be very small. For instance, typical By for the transonic flow simulation is o(10~S- 10-6~. Thus, explicit time integration is not appropri- ate and implicit time integration is indispensable. At the same time, adding high-order artificial dissipa- tion is required to keep the computation stable espe- cially when the central difference is used for convective terms. Second, accurate prediction of aerodynamic co- efficients is required. A few percent of the drag re- duction is important for aircraft design, and the ef- fect of Reynolds number and other parameters should be accurately evaluated. Transition and turbulence modeling are also important. Simulation results for lower Reynolds-number flows can not be extrapolated to high Reynolds numbers. In some applications such as leading-edge separation flow field over a delta wing as shown later, strong separation vortices become key factor for the aerodynamic coefficients. Since the Reynolds number is high, flow field away from the body surface is considered to be rotational inviscid in such applications. However, as is the case of an- other applications, viscous effect near the body surface should be properly evaluated because that is the key factor to determine the location and the strength of the separation vortices. The present author has been engaged in devel- oping an efficient and accurate method for the simu- lation of complicated flows by solving the compress- ible Navier-Stokes equations. In the initial stage of the development, main effort was laid on improving the efficiency. With the progress of supercomputer capability, recent effort tends to be focused on the improvement of the accuracy. In the present paper, the method and the developed Navier-Stokes computer code named LANS3D ~ PLANS' stands for the LU-ADI Navier-Stokes code) is briefly described first. Some of the application examples are then presented. First ex- ample is the transonic flow over transport aircraft con- figuration. Second example is the simulation of vortex breakdown over a strake-delta wing conducted when the author was a research associate at NASA Ames Research Center. In this simulation, the importance of the grid resolution and the accuracy of the numeri- cal algorithm was realized and the computer code was modified. Finally, recent application to the spaceplane configuration is presented. The spacecraft design is so critical that the simulation should be accurate for the s

complicated flow field. These examples will show how the accuracy and efficiency of the code have been im- proved and what problems are left to be improved. 2. Governing Equations and Numerical Algo- r~thm Compressible Navier-Stokes Equations The basic equations under consideration are the unsteady Navier-Stokes equations written for a body- fitted coordinate system ((, 71, () TO + 06E + 0nF + 0<G = Re-l0cS (1) (1) is written as In Eq.~1), the thin-layer approximation has been introduced. The use of thin-layer Navier-Stokes equa- tions is justified because the viscous effects are con- fined to a thin layer near the wall and are dominated by the viscous terms associated with the strain rates normal to the wall, and because the flow away from the body is essentially inviscid. It should be noted that all the viscous terms are not properly evaluated in the viscous layers even by the full Navier-Stokes equations because of the grid deficiency. The pressure, density, and velocity components are related to the energy for an ideal gas by P = flyl)te - 1 pfu2 ~ v2 + wall (2) As will be mentioned later, thin-layer approximation is sometimes extended in two directions. In that case, additional terms appear in Equal) although there is no essential change in the solution process. Numerical Algorithm The algorithm is called LU-ADI (or ADI-LU) factorization method proposed by Obayashi et al.~1] and the 3-D Navier-Stokes code based on this algo- rithm was developed by Fujii and Obayashit2~. Note that this LU-ADI algorithm is different from the orig- inal algorithm presented in Refs. 3 and 4. Since the applications are focused on the Navier-Stokes simula- tions, LU-ADI algorithm uses the Euler implicit inte- gration in time. With the aid of factorization, most of the implicit schemes avoid the memory problem of inverting huge matrices. They are usually classified into two categories based on how to factorize the im- plicit operator. One is the ADI factorization and the other is the LU factorization. In the ADI factoriza- tion method, the left-hand-side implicit operator is split into three components for each direction using approximate factorization. In the LU factorization, the operator is split into two components based on the positive and negative flux Jacobian matrices. The ADI factorization is efficient but has large factoriza- tion error and stability problem, and the LU factor- ization has more arithmetic operations and difficulty for the efficient-vectorization. LU-ADI algorithm tries to decrease the arithmetic operations and can be con- sidered as the compromise of the ADI and LU factor- ization algorithms as pointed out by Jamesont54. In this algorithm, LU factorization is introduced in ad- dition to the ADI factorization,. Each ADI operator which is obtained by the so called Beam-Warming al- gorithmt6] is first diagonalized and then decomposed into two parts using the flux vector splitting ideate. To keep the diagonally dominance, approximate LDU decomposition is used instead of simple LU decompo- . . Sutton. The Beam-Warming factorization applied to Eq. (I + humanDate) (I + h5~Bn _ Dirty x (I + h§<CnhRe-~<J-'MnJDirt) ~Qn =h (Eden + §~Fn + §<Gn _ Re-~<Sn) (Deb + Del,7 + Del<)Qn (3) where h is the time-step, and ~ is a central finite- difference operator. The Di and De terms are implicit and explicit artificial dissipation terms which should be added to the left-hand side and the right-hand side, respectively, to maintain stability. The basic algorithm is first-order accurate in time and second-order accu- rate in space. In the LU-ADI algorithm, each ADI operator is diagonalized and then decomposed into two parts us- ing the technique of flux vector splitting. For example, in the (-direction, I ~ hip A T6(I + ht6D+ + hat DA )T6- · To (IahDAj + h{6D+)(I + ahI DAj 1) x (I + othD+j + half DA )T6-1 where ~ is a coefficient appearing on the j-index for the upwind difFerencing. In other words, cat = 1.0 when first-order differencing is used for the (-derivative, and 1.5 when second-order di~erencing is used. Currently, ~ is set to be 1.0 to eliminate additional arithmetic op- erations. The decomposition in Eq. (4) can be called approximate LDU factorization, or diagonally domi- nant factorization. This decomposition is more stable than simple LU factorization because the diagonal el- ement always has IDA. In the solution process, an inversion in one di- rection consists of one scalar forward sweep, and one scalar backward sweep. Thus, LU-ADI algorithm re- quires little additional memory and is easily vector- ized. Note that an operator in each direction can be considered to be a single iteration of a symmet- ric Gauss Seidel relaxation in one dimension. The DA terms are modified to include an implicit artificial dis- 6

sipation term. Another modification is necessary in the (-direction to evaluate implicit viscous terms. Ad- ditional details of the derivation of the LU-ADI algo- rithm are given in Ref. 1. In the right-hand side, ex- plicit artificial dissipation is necessary. the dissipation model was changed time to time, but currently, non- where linear second-order and fourth-order mixed dissipation model is implemented. The idea comes from the TVD scheme and the model is called simplified TVD-type dissipation. Again, the details are found in Ref. 1. and High-order upwind differencing has become pop- ular in recently-developed Euler methods for com- pressible inviscid flow simulations. This feature has been extended in the straightforward manner for the evaluation of convective terms of Navier-Stokes com- putations. Lately, the matrix form of the dissipation terms implicitly introduced by upwind methods was studded t8,9] and, it was shown that such terms in the upwind schemes such as Roe's flux difference splitting became small in the viscous layers. Thus the use of the proper upwind could improve the accuracy in the viscous layer as well as make the discontinuities sharp. Convective terms in the right-hand side of the code was recently modified and the code includes the option to choose either original central differencing with artificial dissipation or flux difference splitting with Roe's averaget10~. Higher-order extension of flux difference splitting using the MUSCL approach is found in Ref. 8, but is briefly described again. When the convective terms are differenced with the flux-difference splitting of Roe, the spatial derivatives are written in the conservative form as a flux balance. For instance in the (-direction, ~ ~ ~ ABED (Ej+~/2Ej-il2) use The numerical flux Ej+~/2 can be written as the solution to an approximate Riemann problem and the necessary metric terms are evaluated at the cell inter- faces j+1/2. Ej+il2 = 2 tE(Q~ ~ + E(QR)AIRI )]j+~/2' (6) where E is the flux vector and A is the corre- sponding Jacobian matrix computed using the Roe's average state -. QL and QR are the state vari- ables to the left and right of the half-cell interface. These state variables are determined from the locally one-dimensional non-oscillatory interpolations called MUSCL approach. Primitive variables qua, U. V, W. p]T are used for that purpose, and high-order accurate monotone differencing is given by a one-parameter '` . (q~)j+~/2 = qj + 4~1 asps_ + (1 ~ KS)~+]j (qR)j+~/2 = qj+l 4~1 asks+ ~ (1 + ~s)~_]j+l, (~+)j = qj+l - qj, (~_)j = qj - qj_, _ 2~_+6 S (,\+)2 + (a\ )2 + ~ (7) s is the Van Albada's limiter and ~ is a small constant to prevent zero division. For all the results here, third- order accuracy corresponding to arc = 1/3 is used. Near the boundary, the MUSCL interpolation goes back to the first-order. The flux evaluation for the central difference method is written in the same form as Eq.~6) con- sidering the artificial dissipation into account. Ej+il2 = 2 fEj + Ej+, - (~1 - ~52 - ~64V6~+~/2(Qj+i - Qj)] where ~ is a kind of limiter function. The first two terms Ej and Ej+i correspond to the second-order central difference and the last term corresponds to the dissipation model. The same is true for Eq.~6~. The first two terms construct central dif- ference and the last term is dissipative term implicitly added by upwinding. Without limiter functions, first two terms become second-order central difference for first-order upwind method, and become fourth-order central difference for third-order upwind method. The point here is the form of the dissipation terms. In the artificial dissipation model used with the central differ- ence, magnitude of the coefficient is the same for all the equations. Suppose we decompose the Euler equations into a set of independent equations for the waves with the characteristic speed u, u+c, u-c, the magnitude of coefficient should be chosen to be large enough for any of the waves. In the flux evaluation of the Roe's flux difference splitting, this coefficient is in the form of ma- trix A. If you rewrite this dissipation term, it is recog- nized that dissipation depends on the strength of each wave and the magnitude is automatically determined for each wave by upwinding. Thus, artificial dissipa- tion included in the high-resolution upwind methods could be smaller than the dissipation used with the central difference methods. It should also be noted that Vatsa et al. demonstrated this dissipation terms become automatically small in the viscous layers near the body surfaced. On the other hand, some kind of scaling to decrease the dissipation is necessary near the body surface for the central difference approach. 7

3. Application Examples 3.1 Aircraft Simulation - Toward Efficient Code- Following the several trial computations to eval- uate the code capability/4], the developed 3-D Navier- Stokes code based on the LU-ADI algorithm was first used for the simulation of an isolated wing named W- 14~23. This wing geometry was developed as a wing for the transonic transport aircraft by Mitsubishi Heavy Industries. The purpose of this study was to compare the computed result with the experiment, and thus, several computations were carried out for the angles of attack varying from 0° to 7.5°. The Mach number was fixed to be the design Mach number 0.82. Total number of the grid points was 200,000. The angles of attack assigned for the computations were slightly modified from those of the experiment by MHI design- ers, and the computers: Fujii and Obayashi, were not informed of the experimental data in advance. The flow was assumed to be fully turbulent and the so- called Baldwin and Lomax turbulence model[11] was used. The comparison of the surface Cp distribution with the experiment for 2.46- deg. case is shown in Fig. 1. The computed result shows the same tendency as the experiment. The data agreement for both up- per and lower surface is good at every station and, in addition the pressure at the trailing edge is well pre- dicted, which is important for the design process of a new wing. In early 1986, the code was modified for the sim- ulation of wing-fuselage combinationt12~. Again, the used geometry named W-18 was designed by the MHI and seven angle of attacks (from 0.0° to 6.0°) that were shifted a little from the experiment by the empirical way were suggested to the computers who were not informed of the experimental data. Figure 2 shows the grid for the upper half volume of the computa- tional domain. Since not only the viscous layers over the wing but also the viscous layers over the fuselage should be considered, the thin-layer approximation is now adopted for two directions and thus, the basic equations (1) are modified. The turbulence model was also modified near the wing-fuselage junction. The computational grid was generated by Takanashi us- ing his modified conformal mapping technique, and the number of the grid points was increased to about 800,000. Figures 3a-3b show the surface pressure contour plots for the 6.0-de" case. There occurs a large shock- induced separation near the root section, and shock wave exists even at the fuselage surface. The shock wave has a strong spanwise curvature as in the case of isolated wing, but remarkable difference exists near the wing-fuselage junction. For an isolated wing, the shock wave is always perpendicular to the center symmetry plane. On the other hand, the shock wave curves for- ward at the root section because of the large recircu- lating region induced by the wing-fuselage interaction. Computed off-body particle path traces are presented in Fig. 4. The strong spiral is created and moved out- board behind the shock wave over the wing. Near the wing-fuselage junction, there exists a large recirculat- ing region, where the vertical flow resembles the coiled spring bent 90°. The vortex axis is perpendicular to both the wing and the fuselage surface. The Cp dis- tributions over the wing and the fuselage surfaces for 4.0-de" case are compared with the experiment in Figs. 5a and 5b respectively. The overall agreement is fairly good, not only for the wing but also for the fuselage. Note that the discrepancy in the pressure level at the tip section can be explained by the elastic deformation of the test model in the experiment. It is possible that the tip section of the steel model was twisted by aero- dynamic forces, since the aspect ratio of the wing is very high. The discrepancy in the inboard region may be due to the poor turbulence model. The details of the computations are found in Ref. 12. In 1988, along with the parametric study for several wing-fuselage combinations, the flow simula- tion over an almost complete aircraft was tried by Takanashi et alt13~. In addition to the wing and fuse- lage, vertical and horizontal tails were added. The centerline symmetry plane was modified to include the viscous effect of the vertical tail. The horizon- tal tail wing were sandwiched between two chordwise grid lines and the effect of viscous layers over the tail surface was introduced there. Only the coarse grid simulation using about 600,000 grid points was carried out. The surface grid distribution and the computed surface pressure contours are plotted in Fig. 6 and 7, respectively. The Mach number is 0.6, the angle of attack is 0 deg and the Reynolds number is 3.5x106 for this case. Reference 13 shows the details of the computation. All the computations above were carried out within the specified time frame, as they should be for the design purpose. The reliable solutions, thus, had to be obtained within reasonable computer time and the efficiency was the matter of concern. Initially, wing simulation required roughly two hours on Fujitsu VP supercomputer, and now requires 20 min. to 40 min. of computer time. Wing-fuselage simulation could be carried out within 2 to 3 hours. However, for the de- sign purpose, twice as much should be considered for safety. Checking the effect of parameters such as grid spacing or artificial dissipation model is also necessary. 3.2 Vortex Breakdown Simulation - Toward Accurate Code- The flow over aircraft and missiles at moderate 8

to high angles of attack is characterized by the pres- ence of large spiral vortices on the leeward side of the body. These separation vortices induce low pressure on the body upper surface, and this low pressure is the predominant factor of the resulting aerodynamic characteristic of the body. Research on such flow is of great importance practically as well as physically because understanding of the separated and vertical flow fields may lead to the control of vortex behav- ior and eventually to the enhancement of flight vehicle performance. To understand vertical flow field structure over a strake-delta wing configuration, a series of computa- tions was conducted when the author was a NRC re- search associate at NASA Ames Research Centeri14~. Two types of vortex breakdown - bubble and spiral shaped - were successfully simulated and the differ- ence of the flow structure of each breakdown was char- acterized. However, the result at the same time in- dicated still better grid resolution and reducing the artificial dissipation are critically important for an ac- curate simulation of vertical flow field. In all the air- craft simulations shown above, central difference was used for the convective terms. High-resolution upwind difference might have less dissipation as noted in the section 2, and was implemented as an option. Here, some of the results t15] are presented and compared with the result by central difference computation (note that central difference always requires some artificial dissipation model with it). The flow field is the sub- sonic flow over a strake-delta wing. The freestream Mach number is 0.3, and the Reynolds number based on the root chord is 1.3x106 in the following computa- tions. In the first example, the total number of grid points is about 850,000; 119 points in the chordwise direction, 101 point circumferentially and 71 points in the normal direction. Details of the grid genera- tion and the grid distribution can be found in Ref. 14. Figure 8 shows the overall view of the spanwise total pressure contour plots at several chordwise sta- tions at the angle of attack 12 degrees. The contours are plotted at 35% to 95 % chordwise stations with 10 % increase. At this angle of attack, there exist two vortices over the upper surface of the wing; one ema- nating from the strake leading edge and the other from the main-wing leading edge. These two vortices merge together over the main-wing surface because of the mu- tual interaction. Both results indicate the existence of two vortices over the wing surface and their interac- tion. It seems that the merging of the two vortices occurs more downstream in the upwind solution. The corresponding particle path traces showing the vortex trajectories are shown in Fig. 9 for the upwind result. The computed vortex trajectories are presented in Fig. 9. Also presented are the experiment and the result by the central difference computation. It is clear that merging of two vortices occurs downstream in the up- wind result, but still upstream than the experiment at the same Reynolds number. The same computation was carried out using smaller number of grid points (about 120,000 in to- tal). Compared to the previous grid, the number of the grid points are decreased in all the directions. Fig- ures 10a and 10b represent the total pressure contour plots obtained by the upwind and central difference computations, respectively. The contours are again plotted at 35% to 95 % chordwise stations with 10 % increase. The upwind result shown in Fig. 10a indi- cates the existence of two vortices and their merging process although the inner vortex is not as distinct as the fine grid result. On the other hand, the central difference result in Fig. 10b shows only one flattened vortex instead of two vortices. Final example is the result at 30 degrees. The previous study in Ref. 15 showed that vortex break- down takes place near the trailing edge both in the ex- periment and in the computational results on the fine grid. Here medium grid (previously mentioned grid of about 120,000 points) computations are carried out both with the central di~erencing and the upwind dif- ferencing. The computed total pressure contour plots are presented in Figs. Ha and lib. An abrupt in- crease of the vortex-core is observed near the trailing edge in the upwind result shown in Fig. lla. This in- dicates that the vortex has undergone breakdown. In fact, the plot of the streamwise velocity (although not shown here) showed that there exists the reverse flow region near the trailing edge. The central difference result shown in Fig. lib, on the other hand, does not show such a sudden change. Again, the resolution is enhanced by the use of the present upwind scheme at least on the grid used here (although a slight increase of the number of grid points may introduce breakdown phenomenon also in the central difference result). It is recognized from these results that the present upwind scheme has better resolution than the conventional central difference scheme on the same grid although grid resolution itself is, of course, an important factor for an accurate flow simulation. Up- wind scheme is more "vortex-preserving" than central differencing scheme ~ with added dissipation ~ since it has a lower level of dissipation. In the present upwind scheme where the dissipation terms are constructed in the matrix form, each characteristic wave has its own minimum dissipation. On the other hand, cen- tral difference scheme where dissipation terms are con- structed in the scalar form, requires amount of dissi- pation which is large enough for all the waves. Of course, solutions of both central and upwind difference schemes should converge to the solution of 9

the original partial differential equations as the compu- tational grid is refined. The point here is that the accu- racy estimation based on the idea of the Taylor expan- sion is important but not good enough for the system of nonlinear equations. What we need in numerical schemes is the better representation of the properties of original partial differential equations and, in that sense, upwind difference scheme shows better result than that of the central difference scheme for the grid distributions feasible under the memory restriction of the current supercomputers. ~~ ~ ~ A A ~ r 3.3 Spaceplane Simulation - Accuracy and Efficiency Required- Since February 1985, when U.S. President Rea- gan announced the NASP project, there occurred a strong acceleration on the research of orbiting or hy- personic flight vehicle. The CFD is one of the im- portant areas necessary for the development, and the computations have been carried out for transonic to hypersonic flow regime for the spaceplane configura- tion proposed for the research at the NAL in Japan. Following the comparison of the original central di~erencing and the new upwind differencing that is described abovet15], a series of simulations were con- ducted. One of the result is presented in Fig. 12 where the computed surface density contours at the Mach number 1.5, the Reynolds number 4 x106, and the angle of attack 15 ° are plotted. The contours on the wing surface and the fuselage surface indicate the vertical flow generated above. The surface pressure distributions at one cross section is compared with the experiment in Fig. 13. The disagreement in the lower surface Cp at the final station is due to the model sup- port, and otherwise pretty good agreement was ob- tained with the experiment. The details of the series of computations and the effect of each element such as tail fins will appear in Ref. 16. 3.4 Strategy for Complex Configurations One of the important items to be solved in or- der to apply the CFD to practical problems is how to handle complex configurations. Although the so- phisticated grid generation programs have been de- veloped by many researchers, the application of the flow solution codes is still restricted to relatively sim- ple configurations, and a method is needed to make the code applicable to the flow simulation over truly complex configurations (say, complete aircraft with na- celles, engines and so on). Another important item is the problem of the grid resolution. As is seen in the result shown above, even with the high-resolution up- wind scheme, total number of grid points necessary for accurate flow simulations becomes enormous. As the number of the grid points becomes large, computer time also becomes large and, as a result, supercomput- ers with much large memory and faster speed would be required. Recent effort to solve these problems seems to be in two directions. One is the use of unstructured mesh. Using the unstructured mesh, handling complex geom- etry is easier and such grid system is combined with the finite element method or finite volume method for arbitrary cell shapes. So far as the Euler equations are concerned, such approach may be a good choice since explicit time integration can be used. However, for viscous flow simulations, implicit time integration may be necessary and the problem how to apply the implicit time integration and how to adopt turbulence models should be solved. The other approach is the zonal method. In this method, computational grid is constructed for each element of the geometry. For in- stance in aircraft simulation, each of the wing, fuselage and engine has its own grid distributions. Although modification of the existing finite difference codes to the zonal codes is relatively easy, accurate and robust interface method to transfer the data for each zonal region should be developed. In this approach, num- ber of the grid points is also easily increased locally. One of the examplet174 is shown in Fig. 14. The flow field and the solution algorithm are the same as Fig. 10b. Number of the grid points is globally decreased but is locally increase in the vertical flow region, and the result seems to be better than Fig. 10b although the total number of the grid points is almost the same. Currently the effort to develop the zonal interface algo- rithm such as Chimera methodt18] is underway. With such method, the zonal code can handle complex body configuration, and the applicability of the code to the practical problems is to be improved. 4. Supporting System Development With the increase of the data obtained by the CFD research, importance of visualization of the com- puted results began to be recognized. To help people to understand what happens in the flow field from the obtained data, we have to visualize the flow and the use of graphic workstations has an important role for that purpose. Compared to the circumstances sev- eral years ago, the level of the graphic software is im- proved, and visualizing complex flow field in realistic image now is not a difficult task. However, in most cases, such softwares are used only to create nice and beautiful pictures. That may be only important to ad- vertise the CFD capability. What is really important is the development of the software to help understand- ing the result. Such softwares should be interactively used without requiring researcher's effort and can dis- play the plots they want quickly. Displaying beautiful and real-image pictures does not have the first priority. The computer programs that satisfy some of the 0

requirement mentioned above have been developed at NASA Ames Research Center. They are called ' PLOT3D', 'GAS' and 'RIP' and each of them has its own unique featuret19,20~. For instance in RIP, particle trajectories are computed on the CRAY-2 su- percomputer and the result is displayed on the IRIS workstation. The graphic process is carried out inter- actively on the IRIS display. Graphic workstations are also useful for grid gen- eration. Sometimes, grid generation process requires more human time than flow computation itself. Inter- active use of graphic workstation will reduce amount of effort for the grid generation process. One of the super Graphic Workstation has been introduced at our laboratory lately and the preproces- sor (grid generation) and postprocessor (flow visual- ization programs) for Navier-Stokes and other solvers having the features discussed in this paper are in the process of being built up on that machinet214. Here some of the display examples for the post processor are presented. First, it has various function, not only con- tour plotters but also particle path tracer, shock wave detector and so on. Two or more flow functions can be displayed not only on the same window but also on the two or more windows respectively because displaying two or more functions at the same time would help understanding of the physical phenomenon. One of these examples are shown in Figure 15. It shows pres- sure contours of the flow field around a double delta wing (discussed above) on the left window and total pressure contours on the right window. The picture are obtained in the following simple manner. Firstly, the grid and flow data are read in. Then two windows are opened. The left window is selected to display the pressure contours and then the right window is se- lected to display the total pressure contours. The third window can be opened and on which another function can be displayed if necessary. Also the transforma- tions, which are the translation, the rotation and the scaling, can interactively be done as will be described later. Next is the other way of using the multi-windows. We use GWS because output list of enormous numbers obtained by the CFD computation is too difficult to understand. However, after displaying the computed results on the screen we frequently find that we want to extract the digital numbers from the result on the screen. One example is the three-dimensional posi- tions such as the center of a vortex. Figure 16 shows an example of three-dimensional positioning using the multi-windows. The left window is the front view and the right window is the side view of the same objects. Hair cursors move corresponding to the movement of the mouse. If you want to know the position of the vortex center of the third cross section, for instance, it is achieved by setting the hair cursors crossed at the center of vortex on the two windows as shown in Fig- ure 16. Because the right window is perpendicular to the left one, the three-dimensional position is uniquely defined by the two set of the hair cursors. The viewing of any windows can be modified by mouse input while the viewing has its own default val- ues. Users do not need to know the exact definitions of viewing but move the mouse till the desired viewing is achieved. 5. Summary The computer code named LANS3D, one of the representative Navier-Stokes codes in Japan, is taken as an example and the capability of the current CFD technology was discussed. This code was developed for the numerical simulation of high-Reynolds num- ber compressible flows. The algorithm used in this code and how it has been improved so far explained two important aspects of the computational fluid dy- namics (CFD) codes: efficiency and accuracy. Some of the application examples showed the capability of the code for engineering problems. The code and its modified versions have been extensively used by many researchers. Reference 22 reviews such applications. There are many problems (such as turbulence model, unsteady effect) to be solved before making the CFD codes an engineering tool, but imminent problem is the treatment of complex configurations. To use the best of the CFD, supporting system is important and the graphic software development for fluid dynamic re- search on the workstations should have attention. Acknowledgement The LANS3D code was mainly developed by the first author of this paper and Dr. Shigeru Obayashi (currently at NASA Ames Research Center) at the Na- tional Aerospace Laboratory, Japan and many people have helped improving the code capability and its ap- plications. The authors express their special thanks to Dr. Susumu Takanashi and Mr. Masahiro Yoshida at the National Aerospace Laboratory and Ms. Kisa Matsushima at Fujitsu Limited. References 1. Obayashi, S., Matsushima, K, Fujii, K. and Kuwahara, K., 'MImprovements in Efficiency and Reliability for Navier-Stokes Computations Using the LU-ADI Factorization Algorithm, AIAA Paper 86-338", Jan., 1986. 2. Fujii, K. and Obayashi, S., "Navier-Stokes Simu- lations of liansonic Flows over a Practical Wing Configuration," AIAA J., Vol. 25, No. 3, March, 1987 pp. 369-370, also AIAA Paper 86-0513. 3. Obayashi, S. and Kuwahara, K., "An Approx- imate LU Factorization Method for the Com- 11

pressible Navier-Stokes Equations, Journal of Computational Physics, Vol. 63, 1986, pp. 157- 167. 4. Obayashi, S. and Fujii, K., "Computation of Thre - Dimensional Viscous liansonic Flows with the LU Factored Scheme," AIAA Paper 85- 1510, July, 1985. 5. Jameson, A., "Successes and Challenges in Com- putational Aerodynamics," AIAA Paper 87-1184, June, 1987. 6. Pulliam, T. H. and Steger, J. L., "Implicit Fi- nite Difference Simulations of Three-Dimensional Compressible Flow," AIAA Journal, Vol. 18, No. 2, Feb. 1980, pp. 159-167. 7. Steger, J. L. and Warming, R. F., "Flux Vector Splitting of the Inviscid Gas-Dynamic Equations with Application to Finite-Difference Methods," Journal of Computational Physics, Vol. 40, 1981, pp. 263-293. 8. Vatsa, V. N., Thomas, J. L. and Wedan, B. W., "Navier-Stokes Computations of Prolate Spheroids at Angle of Attack," AIAA Paper 87- 2627, August, 1987. 9. Van Leer, B., Thomas, J. L., Roe, P. L. and Newsome, R. W., 'MA Comparison of Numerical Flux Formulas for the Euler and Navier-Stokes Equations," AIAA Paper 87-1104CP, June, 1987. 10. Roe, P. L., "Finite-Volume Methods for the Compressible Navier-Stokes Equations, Proc. Int. Conf. Num. Methods for Laminar and lur- bulent Flows," July, 1987. 11. Baldwin, B. S. and Lomax, H., "Thin Layer Ap- proximation and Algebraic Model for Separated Turbulent Flows," AIAA Paper 78-257, Jan., 1978. 12. Fujii, K and Obayashi, S., "Navier-Stokes Sim- ulations of liansonic Flows over Wing-Fuselage Combination," AIAA J., Vol. 25, No. 12, Dec., 1987, pp. 1587-1596, also AIAA Paper 86-1831. 13. Takanashi, S., Obayashi, S., Matsushima, K. and Fujii, K., "Numerical Simulation of Compressible Flows around Practical Aircraft Configurations," AIAA Paper 87-2410CP, August, 1987. 14. Fujii, K. and Obayashi, S., "Use of High- Resolution Upwind Scheme for Vortical Flow Simulations," AIAA Paper 89-1955, June, 1989. 15. Fujii, K. and Schiff, L. B., "Numerical Sim- ulations of Vortical Flows over a Strake-delta Wing," AIAA Paper 87-1987, 1987, to appear as AIAA Journal in 1989. 16. Matsushima, K.,-Takanashi, S. and Fujii, K., "Navier-Stokes Computations of the Supersonic Flows about a Spac~plane," AIAA Paper 89- 3402, August, 1989. 17. Fujii, K., ' 'A Method to Increase the Accuracy of Vortical Flow Simulations," AIAA Paper 88- 2562CP, 1988. 18. Buning, P. G., Chin, I. T., Obayashi, S., Rizk, Y. M. and Steger, J. L., "Numerical Simulation of the Integrated Space Shuttle Vehicle in Ascent," AIAA Paper 88-4359, 1988. 19. Lasinski, T., et. al., "Flow Visualization of CFD Using Graphics Workstations," AIAA Paper 87- 1180, 1987. 20. Watson, V., et. al., "Use of Computer Graph- ics for Visualization of Flow Fields," presented at AIAA Aerospace Engineering Conference and Show, February 17-19, 1987. 21. Tamura Y. and Fujii, K., "Use of Graphic Work- stations for Computational Fluid Dynamics," to be presented at the International Symposium on Computational Fluid Dynamics-Nagoya, to be held at Nagoya, Aug., 1989. 22. Fujii, K., "Capability of Current Supercomputers for the Computational Fluid Dynamics," to be presented at the 1st International Conference on Applications of Supercomputers in Engineering, to be held at Southampton, U.K. Sept., 1989. 2

°^ EXPERIMENl PRESENT RESULT -1.0 ~ _ ~ L I .o -1 Cp o _ .u 1.0 _ -1.0 -.5 Cp o .5 _ 1.O ) it, ~ I ~6 ~ ' ~ SEMI-SPAN STATION ~ 15.0% _ ~ . ~ _ · · ~ . ~ - ~/-1 ~ ~ 60.0Yo 0 .2 .4 .6 .8 1.0 x/c Wo 37.0% . . . . . _ _. 85.0% . . , . . , .2 .4 .6 .8 1.0 X/C Fig. 1 Comparison of the computed Cp distributions with the experimental data at several spanw~se stations : Moo = 0.82, ct = 2.46°, and Re = 2.0 x 106. Fig. 2 Overall view of the discretized region of the grids (upper half of the computational volume). 13

a) overall view b) close-up view by, "Ye Fig. 3 Computed surface pressure contour plots : Moo = 0.82, ~ = 6.00°, and Re = 1.67 x 106 AT Fig. 4 Computed off-body particle path trace : Moo = 0.82, ~ = 6.00°, and Re = 1.67 x 106. 14

or, deg O EXPERIMENT 3.65 PRESENT RESUl-T 4.00 -1.2 r -.8 -.4 Cp O .4 8 1.2 _ -1.2 -~8 -.4 Cp O .4 .8 1.2 _ ~0Oo _ ,: _ ~ _ SEMI SPAN STATION ~ 13.4% _ ~ .. 1P~ _ ~ _ ~ 60.0X 0 .5 x/c _ - ~~° _ ~ 22.6% ~ , , r Wo ~ ~ '~ 74.6% . . · ~ . 1.0 0 .5 1.0 0 x/c 85.0X .5 1.0 xlc ~ 95.0% . . 0 .5 1.0 x/c Fig. 5a Chordwise Cp distributions over a wing surface at several span-wise stations : Moo = 0.82, ~ = 4.00°, and Re = 1.67 x 106. -1t -.5 ~ ~ . o° P 0f .SI . , , -lr Cp -.5 1 Ol .5 _ 55o o o 1 ~ 0 0 113° ~\ ~f C~o l .35 .40 .45 .50 .55 35 1 39o ~,5~_ . I · . I _ .40 .45 .50 .55 .35 X/L XIL -O ~, 27.57 1 ~ `~55.0 ~_484.32 ~j`~,3 I ~ ' ~ ~139 180 158 G 180 XIL a,de' EXPERIMENT 3.65 ~r .55 Fig. 5b Chordwise Cp distributions over a fuselage surface at several circumferential stations : Moo = 0.82, cx = 4.00°, and Re = 1.67 x 106. 15

Fig. 6 Body geometry and the surface grid for the wing-fuselage-tail combination. Fig. 7 Computed pressure contour plots for the wing-fuselage-tail combination : Moo = 0.60, ~ = 0.0°, and Re = 3.47 x 106. 16

PMAX = 0.9900 P Ml N = 0.8500 :< UP = 0.0050 \~ upwind difference result PMAX = 0.9900 PMIN = 0.8500 \~, UP = 0.0050 \ ~ central difference result Fig. 8 Computed total pressure contour plots at several chordwise stations : Moo = 0.30, ~ = 12.0°, and Re = 1.3 x 106. central difference upwind result result Fig. 9 Computed vortex position compared with experiment : Moo = 0.30, ax = 12.0°, and Re = 1.3 x 106. 17 INNER PRIMARY VORTEX /; OUTER PRIMARY VO RTEX \ EXPERIMENT BY BR ENN ENSTU H L & H UMM E L: Re = 1.3 X 106

PMAX = 0.9900 ~ P Ml N = 0.8300 \\ up = 0.0050 \, PMAX = 0.9900 CUMIN = 0.8500 COP = 0.0050 \~ 0.975 I\ 0.950 , \ I \ '\ upwind difference result central difference result Fig. 10 Computed total pressure contour plots at several chordwise stations -medium grid- : Moo = 0.30, cr = 12.0°, and Re = 1.3 x 106. PMAX = 0.9900 PMIN = 0.7900 UP = 0.010 PMAX = 0.9900 ~,` PMIN = 0.7400 UP = 0.01 ad, <A ~ central difference result upwind difference result Fig. 11 Computed total pressure contour plots at several chordwise stations -medium grid- : Moo = 0.30, ~ = 30.0°, and Re = 1.3 x 106. 18

~- upwind difference result _~ central difference result Fig. 12 Overall view of the computed surface density contour plots over a spaceplane : Moo=15' cx=15°, and Re = 4.0 x 106. i\ ( =\ \ 1 PMAX=~.55\00 ) ) ) ) ~ 1////1 aP0 020 JO/////,/ I i/////// I I upwind difference result Fig. 13 Spanwise surface Cp plots (x/c = 91~o) : Moo=1.5, a=15O, and Re = 4.0 x 106. 19 central difference result

\ :< -A Fig. 14 Computed total pressure contour plots at several chordwise stations-zonal solution- : Moo = 0.30, ~ = 12.0°, and Re = 1.3 x 106. Fig. 15 Pressure and Total pressure contour plots over double delta wing 20

DISCUSSION by F. Stern Most of the effort to improve the calculation results has been devoted to the discretization procedures (de, upwind or central differencing) and grid resolution as opposed to the removal of the thin layer assumption; however, for the complex points considered in the paper, I do not believe that the thin layer approximation is appropriate. Please provide justification and evidence for its case. Also, what modifications were made to the Baldwin-Lomax turbulence model in junction regions? Author's Reply Your question is very important. Please read the paragraph after Eq.(1), P.2 in the manuscript. Computational grids usually used for high Reynolds number flow simulations have very small spacing from the body surface outward but have relatively large spacing along the body surface. In other words, cell aspect ratio is quite large. With such grid, we would never be able to resolve viscous terms associated with the strain rates along the body even though we included all the viscous terms in the program. This was checked by many research scientists at NASA Ames Re- search Center several years ago. The same is true for the region away from the body surface. Since grid spacing is not as fine as that near the body surface, viscous terms are not properly evaluated even though the program computes viscous terms there. Thus, we can say that we are just simulating rotational inviscid flow in the region away from the body surface even when we are solving the Navier-Stokes equations. Again, I have checked how large the contribution by the viscous terms is compared to other terms, and the result said no contribution. However, because of the artificial dissipation, the solution tends to be dissipative there. For instance, strength of the vortex tends to decay easily even with the Euler computations. Thus before trying to evaluate viscous terms properly, we have to minimize the artificial dissipation effect and this can be realized by the use of high-resolution upwind method. Inclusion of all the viscous terms is not a difficult task and probably would require only 10 to 20% more computational time. Number of computational grid points is becoming larger and larger. With such background, I think it is the time to evaluate all the viscous terms to remove ambiguity. Thin-layer approximation is only a practical approximation and should be given up soon with the increase of the computational grid points. Please refer [A1] about modification of turbulence model at juncture. [A1] Hung, C.M. and Buning, P.G., "Simulation of Blunt-Fin Induced Shock Wave and Turbulent Boundary Layer Interaction", J. Fluid Mech., Vol.154, pp.163-185, 1985. DISCUSSION by Y. Kodama You showed comparison of "coarse grid" case and "fine grid" case. But in some cases "fine grid" is not fine enough to resolve the flow field. When you compute a flow, on what principle do you determine the number of grid points? Author's Reply The best way to determine the number of the grid points is to compute the flow field until the solution does not change with further increase of the grid points. However, in practice, it is impossible to do so especially in the case of large scale com- putation. The fine grid I used in the computation was the maximum number of the grid points that fitted into the main memory of the super-computer used for that computation. We considered that the solution resolved at least the global feature of the flow by comparing it with the experiment. It is true that the solution I have shown lacks capturing small scales and thus, did not go into the quantitative comparison of the detail (See AIAA J. Sept., 1989 by Fujii and Schiff). Besides, there were many more assumptions such as flow symmetry (which never realize when breakdown takes place) and no turbulence model. We can still say that vertical flow requires fine grid resolution from the fact that the solution depended on the grid resolution. Further research in the near future using more grid points will show how reasonable the current solution is. 22