Solution of aerospace problems using structured and unstructured strategies

Solution of Aerospace Problems using Structured and Unstructured Strategies

Edisson Sávio de Góes Maciel

Instituto Tecnológico de Aeronáutica

Centro Técnico Aeorespacial

12228-900 São José dos Campos, SP. Brazil

[e-mail]

João Luiz Filgueiras de Azevedov

Instituto de Aeronáutica e Espaço

Centro Técnico Aeorespacial

12228-904 São José dos Campos, SP. Brazil

[e-mail]

Products developed at industries, institutes and research centers are expected to have high level of quality and performance, having a minimum waste, which require efficient and robust tools to numerically simulate stringent project conditions with great reliability. In this context, Computational Fluid Dynamics (CFD) plays an important role and the present work shows two numerical algorithms that are used in the CFD community to solve the Euler and Navier-Stokes equations applied to typical aerospace and aeronautical problems. Particularly, unstructured discretization of the spatial domain has gained special attention by the international community due to its ease in discretizing complex spatial domains. This work has the main objective of illustrating some advantages and disadvantages of numerical algorithms using structured and unstructured spatial discretization of the flow governing equations. Numerical methods include a finite volume formulation and the Euler and Navier-Stokes equations are applied to solve a transonic nozzle problem, a low supersonic airfoil problem and a hypersonic inlet problem. In a structured context, these problems are solved using MacCormack’s implicit algorithm with Steger and Warming’s flux vector splitting technique, while, in an unstructured context, Jameson and Mavriplis’ explicit algorithm is used. Convergence acceleration is obtained using a spatially variable time stepping procedure.

Keywords: Euler and Navier-Stokes Equations; nozzle, airfoil and inlet llows; MacCormack scheme; Jameson and Mavriplis scheme; flux vector splitting

Introduction

In the aerospace and aeronautical industries, the need for practical tests in several aerodynamic components of airplane and aerospace vehicle during the project phase is restricted by the high costs to manufacture scaled models and to perform the wind tunnels tests. Other main difficulties are related to the large number of experimental tests required during optimization of these models. The development of computer technology, allowing the existence of machines with high speed processors and high storage capacity, has boomed CFD towards a meaningful role in several sectors of the industry. Such sectors require low level of experimental development costs and products yield with the desired performance.

In this context, this work presents some numerical algorithms used in CFD community to solve the Euler and Navier-Stokes equations applied to traditional aerospace and aeronautical problems. Transonic nozzle, low supersonic airfoil and hypersonic inlet problems are studied in the context of a structured and an unstructured discretization of the flow governing equations. The latter is usually recommended instead of the former, for complex configurations, due to the ease and efficiency with which such domains can be discretized (Jameson and Mavriplis, 1986, Mavriplis, 1990, Batina, 1990, and Pirzadeh, 1991). However, the issue of unstructured mesh generation will not be discussed in this work. Techniques such as Delaunay triangulation (Hefazi and Chen, 1992, and Mavriplis, 1995), which represents an unique triangulation of a given set of points that exhibits a large class of well-defined properties, and the advancing front method (Pirzadeh, 1991, Mavriplis, 1995, and Korzenowski and Azevedo, 1996), which begins with a discretization of the geometry boundaries as a set of edges and advances out these ones into the field, are currently used in the CFD community. The main objective of this work is to highlight the numerical features of MacCormack’s (1985) and Jameson and Mavriplis’ (1986) algorithms in the solution of the Euler and Navier-Stokes equations, regardless of the method used for grid generation.

In the present paper, MacCormack’s implicit scheme (1985) with Steger and Warming’s flux vector splitting (1981) and Jameson and Mavriplis’ explicit scheme (1986) are compared, both in a cell centered type discretization, using a finite volume formulation. The MacCormack scheme (1985) was an excellent option to obtain high convergence rates in the solution of initial value problems. The Jameson and Mavriplis scheme (1986) is the most currently employed scheme in terms of unstructured discretization of flow governing equations (Mavriplis and Jameson, 1987, Batina, 1990, Arnone, Liou and Povinelli, 1991, Long, Khan and Sharp, 1991, Swanson and Radespiel, 1991, and Hooker, Batina and Williams, 1992). The present comparison intends to emphasize important features of these numerical schemes in the following topics: computational performance, some aspects of solution quality and robustness properties.

Nomenclature

Latin Letters

a_* X= critical speed of sound, m/s

a = speed of sound in fluid, m/s

CFL = “Courant-Friedrichs-Lewy” number

E_e = inviscid flux vector (or Euler flux vector) in x direction

E_v = viscous flux vector in x direction

e_i = specific internal energy of fluid, J/kg

e = total energy of fluid per unity volume, J/m3

F_e = inviscid flux vector (or Euler flux vector) in y direction

F_v = viscous flux vector in y direction

P_t = stagnation pressure, N/m²

p = static pressure of fluid, N/m2
q = velocity vector intensity, m/s
u = x component of velocity vector q, m/s
v = y component of velocity vector q, m/s

Greek Letters

a = angle of attack, degrees

g = ratio of specific heats, adopted 1.4 to atmospheric mean

l = diagonal matrix of eigenvalues associated to Euler equations in two-dimensions, according to MacCormack (1985) interpretation of normalised area vectors

r = fluid density, kg/m3

t = viscous stress, N/m²

Subscripts

e = Euler

i = internal or computational index i

j = computational index j

v = viscous

x = spatial position x in Cartesian coordinate system, m

y = spatial position y in Cartesian coordinate system, m

_h = curvilinear direction normal to geometry

_x = curvilinear direction tangent to geometry

Superscripts

ⁿ = counter of time iterations

MacCormack Algorithm

Spatial and Temporal Discretization

The Navier-Stokes equations, in conservative integral form and in a generalized curvilinear coordinate system, can be written as

being the vector of conservative variables, constant in each mesh cell. Mathematically, it is written as

where V is the area of the computational cell (two-dimensional case); represents the gradient operator; and, represents the Navier-Stokes flux:

with:

In the above equations, r is the fluid density; u and v are Cartesian components of the velocity vector in the x and y directions, respectively; e is the total energy per unit fluid volume; p is the static pressure; t's are the viscous stresses; m is the fluid dynamic viscosity; g is the ratio of specific heats; and Prd is the Prandtl number. The Reynolds number characterizing the flow regime is represented by:

where l is a characteristic length of the geometry under study and u_ref is the reference speed.

As the Euler equations are simply obtained by neglecting viscous terms in Eq. (1), the following theory will be presented for the general case of Navier-Stokes equations. In these equations, the derivatives present in viscous terms are calculated by applying Green’s theorem. In other words, such derivatives are considered constant for each volume and they are calculated as a surface integral (Mavriplis and Jameson, 1990, and Long, Khan and Sharp, 1991, and Kwon and Hah, 1995).

The experience has shown that different nondimensionalization are better suited for the different physical problems here analyzed. Hence, for the case of transonic nozzle, Eq. (1) is nondimensionalized with respect to stagnation properties whereas for the case of low supersonic airfoil and for hypersonic inlet in relation to freestream properties. In each case, density is nondimensionalized in terms of stagnation or freestream density; u and v velocity components are nondimensionalized in relation to the speed of sound; pressure and total energy per unit of fluid volume are nondimensionalized in relation to the product between stagnation or freestream density and speed of sound squared. Other details of the present nondimensionalization may be found in MacComack (1984) and Jameson and Mavriplis (1986).

Applying Green’s theorem to Eq. (1) and adopting a structured mesh notation, it is written that:

where n is the outward normal vector for each flux surface with unit norm. In this notation, each computational cell is defined by nodes (i,j), (i+1,j), (i+1,j+1) and (i,j+1). The time integration is accomplished using the following first order scheme:

where the overbar was omitted in the vector of conservative variables for clarity. The discretization of the surface integral in Eq. (6) results in:

where, for example, S_i,j–1/2 has the same direction and orientation as n_i,j–1/2 and norm equal to S_i,j–1/2.

The subscripts indicate that fluxes are calculated in each respective cell surface. The flux areas are determined by S = S_xi + S_yj and are calculated according to Table 1.

The cell areas are determined by the expression below:

In a generalized curvilinear coordinate system, MacCormack (1984) suggests the use of normalized area vectors, orientated in the positive x and h directions. These vectors are defined as:

where . The s_x and s_y terms, associated with the Sx and Sy terms presented above, are defined in Table 2.

A good advantage of this procedure is the computational economy in area vector calculations. It is only a matter of calculating vectors in faces (i+1/2,j) and (i,j+1/2), or (i-1/2,j) and (i,j-1/2), and add or subtract an unit in appropriating subscripts to determine the desired flux term. Equation (7) can be rewritten using the following expressions

The MacCormack schemes are Lax-Wendroff type methods rewritten in the predictor/corrector context, in which space and time are discretized in conjunction with the same order of accuracy. In this procedure, a certain amount of dissipation appears implicitly in the resulting algorithm. Literature commonly references MacCormack schemes as upwind. As a matter of fact, it is only the MacCormack scheme (1985) using flux vector splitting, which includes an analysis of propagation of information in characteristic directions, that can be classified as an upwind scheme. The MacCormack scheme (1985) takes into account Steger and Warming’s flux vector splitting (1981). In accordance with MacCormack’s methodology, convective fluxes in face (i+1/2,j) are written as:

For example, the convective fluxes are calculated as

By definition, the reconstructed Jacobian matrices can be determined as follows:

The same reasoning is true to other flux faces. In the Steger and Warming’s flux vector splitting (1981), it is verified that these positive and negative reconstructed matrices can be determined as:

where matrices T_x , , T_h and are defined, for example, in Chaussee and Pulliam (1981). In this case, based on the equivalence between finite difference and finite volume formulations, in generalized curvilinear coordinates, it is possible to obtain for the predictor step

· Matrices

and :

and , where and are parameters to construct matrices T_x and . Eigenvalues of Euler equations in x direction are:

and

· Matrices

and :

(T_h and T_h^-1) ® k_x = , k_y = , _x = , _y = and where k_x,k_y,_x and _y are parameters to construct matrices T_h and . Eigenvalues of Euler equations in h direction are:

and

where "a" is the speed of sound. Calculating the convective flux terms of Eq. (7) according to Eqs. (10) through (14) and determining the reconstructed Jacobian matrices like (15) and (16), the MacCormack scheme (1985), in its explicit version, is written as:

· Predictor step:

· Corrector step:

The MacCormack implicit scheme (1985) uses the development described above with flux vector splitting to obtain a five-diagonal full block linear system. This system presents a dominant main diagonal when compared to the others. With this procedure, when large time steps are considered, the resulting algorithm is closer to the Newton iterative method and an excellent convergence rate to the steady state solution will be reached. MacCormack’s implicit scheme (1985) is obtained using the implicit Euler method to accomplish time integration in both predictor and corrector steps, resulting in:

· Predictor step:

· Corrector step:

It is noted that the above five-diagonal system has a dominant main diagonal when compared to the adjacent diagonals (MacCormack, 1984 and 1985). This is the main advantage of MacCormack’s implicit scheme (1985). One disadvantage is that, with this implementation, MacCormack method (1985) becomes a first-order scheme in space.

This five-diagonal system, in both predictor and corrector steps, is solved by a symmetric-line-Gauss-Seidel relaxation method with two sweeps: one forward and other backward. This procedure results in a three-diagonal linear system in the h direction. The first time step adopts the trivial solution for the correction in the conserved variables when solving the predictor step. In the corrector step the currently available solution for correction is used. Afterwards, all other Gauss-Seidel sweeps use the latest available values for the corrections. According to MacCormack (1984), the symmetric-line-Gauss-Seidel relaxation method should be solved only with two sweeps to avoid excessive computational cost.

For the sake of computational cost reduction and taking into account that the present work is interested in steady state solutions, the viscous Jacobian matrices at the left-hand side of Eqs. (19) and (20) were neglected. This procedure was also adopted by Chaussee and Pulliam (1981), due to the same reasons, and no deterioration in terms of quality of their results was observed.

Jameson and Mavriplis Algorithm

Spatial and Temporal Discretization

The Navier-Stokes equations in conservative integral form and in a finite volume formulation can be written, in an unstructured context and after spatial discretization, as:

where is the discrete approximation of flux integral in Eq. (21). In the present work, it is assumed that:

with i indexing a given mesh volume and k its respective neighboring volume; n1 and n2 represent consecutive nodes of the i-th volume in counter-clockwise orientation.

The spatial discretization proposed by Jameson and Mavriplis (1986) is equivalent to a second-order centered scheme in a finite difference context. In this way, it is necessary to explicitly introduce an artificial dissipation operator "D" to avoid, for example, odd-even uncoupled solutions and nonlinear instabilities (shock waves). Then, Eq. (21) is rewritten as

The time integration is accomplished by using a second-order, explicit, Runge-Kutta method with five stages that can be generally represented by:

with k = 1,..., 5; m = 0, 2 and 4; a 1 = 1/4, a 2 = 1/6, a 3 = 3/8, a 4 = 1/2 and a 5 = 1. According to Swanson and Radespiel (1991), the artificial dissipation should only be evaluated in odd stages, aiming CPU time economy and better stability conditions based on the hyperbolic/parabolic features of the Navier-Stokes equations. For the Euler equations, m = 0 (k = 1) and 1 (k = 2). The dissipation operator is "frozen" for the reminiscent stages, exploring the hyperbolic features of these equations to assure stable convergence.

Artificial Dissipation Operator

The artificial dissipation operator presented in this work is suggested by Mavriplis (1990) and it has the following structure:

where:

, named undivided Laplacian operator, gives numerical stability in presence of shock waves; and , named biharmonic operator, is responsible for background stability (for example, odd-even uncoupled instabilities). In this last term, . Whenever k represents a special boundary cell, named "ghost" cell, its contribution to the summatory is neglected. The e terms are defined as follows:

with representing a pressure sensor to identify regions of high gradients. The K(2) and K(4) constants have typical values of 1/4 and 3/256, respectively. Again, whenever k represents a ghost cell, it is assumed n_g = n_i. The Ai terms are contributions of the maximum normal eigenvalue of the Euler equations integrated along each cell face. These terms are defined as:

where u_i,k, v_i,k and a_i,k are calculated as the arithmetical average between the respective property values associated with the i-th real volume and its k-th neighboring volume.

Spatially Variable Time Step

With the purpose of accelerating MacCormack’s implicit scheme (1985) and Jameson and Mavriplis’ explicit scheme (1986), a spatially variable time step in each computational mesh cell is used. The basic idea of this procedure is to maintain a constant CFL number in the overall calculation domain, allowing the use of appropriating time steps for each specific mesh region during the convergence process. In this way, according to CFL definition, it is possible to write:

where CFL is the "Courant-Friedrichs-Lewy" number to provide numerical stability; (Ds)_cel: is a characteristic length of information transport. In a finite volume formulation, (Ds)_cel: is chosen as the smallest among the smallest cell centroid-neighboring centroid distance and the smallest cell side length. In the above equation, is the maximum characteristic speed of information transport.

Initial Condition

The initial condition values adopted for the transonic nozzle problem are those of stagnation properties applied to the overall calculation domain, except at the exit boundary where pressure and density are reduced by a given factor to initialize the flow simulation (MacCormack, 1985, Maciel and Azevedo, 1997, and Maciel and Azevedo, 1998). The flow entrance angle is 0.0° and the vector of conservative variables is defined as

· Overall domain except at the exit boundary:

· Exit boundary:

For the other physical problems, freestream values are adopted for all properties as initial condition in the overall domain (Jameson e Mavriplis, 1986). Therefore, the vector of conservative variables is

where M_¥ is the freestream Mach number, a is the flow angle of attack and g is the ratio of specific heats.

Boundary Conditions

The boundary conditions are basically of four types: wall, entrance, exit and symmetry. These boundary conditions are implemented, as commented before, in ghost cells:

a) Wall condition: For the Euler equations, wall condition implies in flow tangency. It is accomplished by considering the velocity component of the ghost volume tangent to the wall to be equal to the corresponding velocity component of its neighboring real volume. At the same time, the velocity component of the ghost volume normal to the wall is taken to be equal in value but of opposite sign relative to the velocity component of its neighboring real volume.

The fluid pressure gradient normal to the wall is assumed equal to zero in the Euler formulation. The same hypothesis is applied to the temperature gradient. From these assumptions, ghost volume pressure and density are extrapolated from its neighboring real volume (zero-order extrapolation).

The non permeability and adherence conditions are imposed for flow at the wall for the Navier-Stokes equations. This is done by imposing that ug = -ur and vg = -vr, where "g" represents the ghost volume and "r" represents its neighboring real volume. The same conditions, as in the Euler case, are imposed for pressure and temperature gradients normal to the wall, according to the boundary layer theory and adiabatic wall condition for the viscous case, respectively.

b) Entrance Condition:

b.1) Subsonic Flow: Three property values need to be specified at this boundary, based on an analysis of propagation of information in the characteristic directions in the calculation domain (Azevedo, 1992, and Maciel and Azevedo, 1998). In other words, for subsonic flow, three characteristic lines have direction and orientation pointing inward to the calculation domain and should be fixed. Only the characteristic line of speed "(qn-a)" cannot be fixed and should be determined by interior information. For the nozzle problem, the velocity component "u" of the ghost volume is extrapolated from its neighbor. In the other problems, pressure of ghost volume is extrapolated from its neighbor.

b.2) Supersonic Flow: all variables are fixed with freestream values at the entrance boundary.

c) Exit Condition:

c.1) Subsonic Flow: Three characteristic lines have direction and orientation pointing outward from calculation domain. These ones are extrapolated from interior domain and the characteristic line of speed "(qn-a)" should be specified. Then, pressure of the ghost volume is fixed at its initial value.

c.2) Supersonic Flow: In this case, all variables are extrapolated from interior domain because the four characteristic lines of the Euler equations are pointing outward from interior domain and nothing can be fixed at the boundary.

d) Symmetry Condition: This condition only exists for the nozzle problem and ghost cells should have a vector of conservative variables defined as below:

because the mesh only discretizes half nozzle. This condition guarantees no crossflow at center line.

Results

The computational performance, general aspects of solution quality and robustness characteristics of MacCormack’s implicit scheme (1985) and Jameson and Mavriplis’ explicit scheme (1986) are presented in the forthcoming discussion. Tests were performed in a PENTIUM-200MHz microcomputer, 64 Mbytes of RAM and the criterion adopted to obtain a converged solution was the order of magnitude of the maximum residue of all conservation equations to be less than or equal to 10^-6.

In the viscous case, a turbulence model was not implemented and, therefore, all tests considered laminar flow. This hypothesis is not a precise description of fluid flow due to the high Reynolds numbers studied herein. However, the authors’ intention is to gain experience in treating flows governed by Navier-Stokes equations in laminar regime. Turbulent flow simulations are beyond the scope of this paper and they represent a future step in the group’s development process.

Unstructured meshes were created transforming each rectangular cell of given structured meshes into two triangular cells. All necessary tables were generated and a volume-based data structure was implemented. Although this procedure of mesh generation does not produce meshes with the best spatial discretization, meshes with reasonable quality have been obtained for the present problems.

Transonic Nozzle

Euler Equations

An algebraic mesh with 760 rectangular real volumes and 820 nodes (41-x x 20-h ) was used for the transonic nozzle problem. An exponential stretching of 10% near the nozzle wall and near the throat was used. The angle of inclination of the nozzle’s wall convergent region is 22.33° and the angle of inclination of the divergent wall region is 1.21° . The nozzle’s total length and its curvature ratio are 0.1158m (0.380ft) and 0.0274m (0.090ft), respectively. The throat is located at the nozzle’s center and its half height is 0.0137m (0.045ft). Figures 1 and 2 show the structured and unstructured meshes generated for the Euler solutions. Figures 3 and 4 show pressure and Mach number contours obtained by MacCormack’s implicit scheme (1985). A CFL number of 125.0 was used. The Computational cost per iteration and per cell was 0.00150s.

Figures 5 and 6 show the same type of solutions obtained by Jameson and Mavriplis’ explicit scheme (1986). The maximum CFL number of 0.9 was used and its computational cost per iteration and per cell was 0.00015s. Comparing similar solutions, it is noted that MacCormack’s scheme (1985) cannot detect the weak shock wave in pressure contours (Fig. 3) near the throat, whereas the Jameson and Mavriplis’ scheme (1986) can efficiently detect it. The excess of numerical diffusion originated by the first-order accuracy in space of MacCormack’s scheme (1985) is the main reason for smoothing the weak shock wave in the pressure contours. Such behaviour of generating excess of diffusion and loss of accuracy in the solution, associated with schemes of first-order spatial accuracy, is well known from literature: Harten (1983), Harten (1984), Sweby (1984) and Liou (1996). In spite of the first-order accuracy of MacCormack’s scheme (1985), Mach number contours agree very well as compared as Jameson and Mavriplis’ solution.

Figures 7 and 8 show convergence histories of each method obtained for this problem. The convergence rate of MacCormack’s scheme is excellent. The steady state solution is obtained in less than 100 iterations due to the appropriated consistency between left- and right-hand sides.

Navier-Stokes Equations

An algebraic mesh of 2,340 rectangular real volumes and 2,440 nodes (61-x x 40-h ) was used for viscous case. The same exponential stretching as for the inviscid case was implemented. The Reynolds number was estimated from the initial conditions, using values of the stagnation properties. Values of density, speed of sound and viscosity were determined according to Fox and McDonald (1988) and assuming a temperature of 294.8 K. The reference speed in this case is taken as the speed of sound at this temperature, resulting in , where g assumed the value 1.4. Hence, the Reynolds number based on a reference length of 0.123 ft is

.

The Prandtl number (Prd) was assumed as 0.72 and the molecular viscosity (m) was assumed constant and equal to 1.0 by nondimensionalization. Meshes are shown in Figures 9 and 10 for both structured and unstructured contexts.

Figures 11 and 12 show pressure and Mach number contours obtained by MacCormack’s method (1985) and Figures 13 and 14 show the respective contours obtained by Jameson and Mavriplis’ method (1986). A CFL number of 20.0 was used by the MacCormack scheme (1985) because of the lack of consistency between right- and left-hand side operators in the viscous case. However, this scheme still maintains a better computational performance than the Jameson and Mavriplis scheme (1986). A CFL number of 1.1 was used by the latter and the steady state was reached in more than 20,000 iterations while MacCormack’s scheme (1985) needed less than 4,000 iterations to obtain convergence. The MacCormack scheme (1985) cannot detect the weak shock wave near the throat, even with a more refined mesh. Solutions are in good agreement in terms of Mach number contours.

Figure 15 shows pressure ratio distribution, in relation to stagnation pressure P_t, along the nozzle wall. The weak shock wave is well predicted at the throat by both schemes, emphasizing the good approximation of MacCormack’s method (1985). Near the nozzle entrance, the results obtained with MacCormack’s scheme (1985) are worse than those obtained with Jameson and Mavriplis’ scheme (1986), evidencing a loss of precision in this region. It is again because of its first-order precision in space. It should be emphasized that quantitative and qualitative results of Jameson and Mavriplis’ scheme (1986) are in good agreement with results presented in Maciel and Azevedo (1997) and Maciel and Azevedo (1998). The computational cost per iteration and per cell for MacCormack’s scheme (1985) was 0.00170s while it was 0.00056s for the Jameson and Mavriplis’ scheme (1986) in the viscous case. It is noted that, while MacCormack’s scheme (1985) is 10 times more expensive than Jameson and Mavriplis’ scheme (1986) in the inviscid case, this difference is only about 3 times in the viscous case. In Jameson and Mavriplis’ scheme (1986), gradients of primitive variables are calculated in the five stages of the Runge-Kutta method. In MacCormack’s scheme (1985), such calculations are only accomplished in two stages (predictor and corrector steps). This is the main reason for such differences in terms of computational cost.

Figures 16 and 17 show convergence histories. The best stability and consistency of MacCormack’s implicit scheme (1985) is again highlighted. MacCormack’s scheme (1985) is about 6 times faster than Jameson and Mavriplis’ scheme (1986).

It is possible to conclude in this problem that MacCormack’s scheme (1985) is an excellent numerical tool to obtain fast results. This type of characteristic is important in industrial applications.

Low Supersonic Airfoil

Euler Equations

An "O" type mesh with 2,352 rectangular real volumes and 2,450 nodes (49-x x 50-h ) was generated for the airfoil configuration. Unit dimension was adopted for the airfoil chord. The far field was located at 25 chords in relation to airfoil leading edge. Initial conditions were set as freestream Mach number equals to 1.0 and the angle of attack equalled to 0.0° for a NACA 0012 profile.

Figures 18 and 19 show structured and unstructured meshes used for this physical problem. Figures 20 and 21 show pressure and Mach number contours obtained by MacCormack’s scheme (1985). A CFL number of 100.0 was used. A fast convergence rate was obtained, although it was worse than that obtained by the nozzle problem.

Figures 22 and 23 show the results obtained by Jameson and Mavriplis’ scheme (1986). A CFL number of 0.9 was employed. It is noted, in this second example, a larger difference, both qualitatively and quantitatively, between the results obtained with the two different schemes. Jameson and Mavriplis’ scheme (1986) underpredicts the pressure distribution and overpredicts the Mach distribution in relation to MacCormack’s scheme (1985). Results obtained with the Jameson and Mavriplis’ scheme (1986) are considered more precise than those obtained with MacCormack’s scheme (1985) because of the second-order of accuracy in space of the former. In terms of pressure, it originates a significant discrepancy on the calculation of lift and moment aerodynamic coefficients. Jameson and Mavriplis’ scheme (1986) detects a greater Mach number peak than MacCormack’s scheme (1985) near the trailing edge. Again, the first-order accuracy of MacCormack’s scheme seems to be the reason for such behaviour.

Figures 24 and 25 show convergence histories. MacCormack’s scheme (1985) reached a steady state solution in less than 200 iterations while Jameson and Mavriplis’ scheme (1986) reached the same order of magnitude of the residue in approximately 16,000 iterations. Jameson and Mavriplis’ scheme (1986) presents more oscillations in the maximum residue history during the convergence process. One reason for this problem is mainly due to the high aspect ratio of triangles near the trailing edge. This high deformation of triangles originates large differences in volume calculation (or Jacobian calculation) near the airfoil which represent severe changes associated with spatial coordinate transformation. In other words, the flux balance at triangles near the trailing edge will not be satisfied so easily during the time march. As consequence, the trailing edge region takes longer to converge to steady state than the rest of the flowfield. This problem is more severe for viscous cases due to the need of capturing boundary layer effects.

Navier-Stokes Equations

An "O" type mesh with 2,832 rectangular real volumes and 2,940 nodes (49-x x 60-h ) was generated. A stretching factor of 10% was used in h direction on the airfoil surface. The far field was also located at 25 chords. The Reynolds number was estimated based on a flight altitude of 10,000m and considering an airfoil chord of 1.0 m as the reference length. The reference speed is the local speed of sound at this altitude, which yields u_ref = c = 299.5 m. Hence, the Reynolds number is

The same values of Prd and m were used as in the transonic nozzle problem. Figures 26 and 27 show structured and unstructured meshes used for the viscous case.

Figures 28 and 29 show pressure and Mach number contours obtained by MacCormack’s scheme (1985), while Figures 30 and 31 correspond to the results obtained with Jameson and Mavriplis’ scheme (1986). CFL numbers of 7.0 and 0.3 were used with MacCormack’s and Jameson and Mavriplis’ schemes, respectively.

The K⁽²⁾ and K⁽⁴⁾ dissipation coefficients were modified in this problem in relation to their standard values for the Jameson and Mavriplis’ scheme (1986) because this problem is more severe than the transonic nozzle problem. Values of 0.4 and 0.03 were adopted, respectively. This amount of artificial dissipation only aimed to guarantee a stable convergence without degenerating the solution quality.

In the above figures, some differences are noted in terms of solution quality. Jameson and Mavriplis’ solution is non symmetrical in relation to airfoil chord direction because of the lack of symmetry of the mesh. Triangles are constructed oriented in the same direction in the overall mesh. In quantitative terms, the value of pressure peak at leading edge in Jameson and Mavriplis’ scheme (1986) is reduced as compared as MacCormack’s scheme (1985). Again, this type of behaviour, also verified in inviscid case, originates errors in the calculation of aerodynamic coefficients. An adequate second-order implementation of MacCormack scheme (1986) can resolve this problem. A suggestion to obtain such accuracy is the use of a "MUSCL" extrapolation procedure (van Leer, 1979).

Figures 32 and 33 show convergence histories in the viscous case. In less than 800 iterations, the steady state solution was obtained by MacCormack’s scheme (1985), while it takes more than 50,000 iterations for Jameson and Mavriplis’ scheme (1986) to satisfy the same convergence criterion. MacCormack’s scheme (1985) is about 60 times faster than Jameson and Mavriplis’ scheme (1986) in terms of number of iterations to reach steady state. Fast convergence is the main important feature of MacCormack’s scheme (1985) and, although the scheme has only first-order accuracy in space, its results describe appropriately the main characteristics of fluid flow. The excellent convergence rates obtained for these two physical problems, transonic nozzle and low supersonic airfoil, highlight more advantages in using this scheme than Jameson and Mavriplis’ scheme (1986).

Particularly, significant oscillations in Jameson and Mavriplis scheme (1986) was observed when the residue is reduced about 4 orders of magnitude. The above convergence history was finished after 50,000 iterations. The results indicate that steady state is reached after 20,000 iterations because no meaningful variations occur in terms of solution features. The reason for this instability is due to the great triangle’s geometric deformation, near the airfoil, as mentioned in the inviscid case.

Hypersonic Inlet

Euler Equations

An algebraic mesh with 2,691 rectangular real volumes and 2,800 nodes (70-x x 40-h ) was created. The initial condition was set as freestream Mach number of 15.0 and angle of attack of 0.0° . Figures 34 and 35 show structured and unstructured meshes used.

Figures 36 and 37 show pressure and Mach number contours obtained by MacCormack’s scheme (1985). A CFL number of 10.0 was employed. Figures 38 and 39 show results generated by Jameson and Mavriplis’ scheme (1986) using a CFL number of 0.9. In qualitative terms, there are good agreement between the results. However, in quantitative terms, there are some discrepancies. MacCormack’s scheme (1985) underpredicts the pressure field, opposed to the behaviour identified in the last problems, while Jameson and Mavriplis’ scheme (1986) predicts a pressure field more severe. This behaviour of MacCormack’s scheme (1985) can again be attributed to its excessive diffusive characteristics due to the first-order spatial accuracy. It is mainly noted in problems with severe shock waves, where more attenuation of physical instabilities is required. Iterations of lower and upper shock waves were well captured. Expansion and compression regions are appropriately described and low level of oscillations in the Jameson and Mavriplis scheme (1986), a second-order method, are obtained mainly due to the features of the dissipation operator. Mavriplis’ artificial dissipation operator (1990) presents good behaviour allowing a more adequate way to preserve propagation of information and to solve instabilities deriving from errors of discretization. Standard values were used for the dissipation coefficients in Jameson and Mavriplis’ scheme (1986) without generating difficulties to obtain convergence.

Figures 40 and 41 show convergence histories. MacCormack’s scheme (1985) reached steady state in less than 300 iterations. Jameson and Mavriplis’ scheme (1986) needed more than 20,000 iterations. It is important to note in this result that steady state solution was obtained by MacCormack scheme (1985) in approximately 200 iterations, where some numerical instabilities occurred, because of the maximum value of CFL number used, and originated deterioration in the convergence rate. In practical studies, a reduction in four orders of magnitude of the maximum residue is more than satisfactory to guarantee a converged solution. Note that, for this hypersonic problem, practical good results are provided by MacCormack’s scheme (1985) in about 150 iterations. As a conclusion, this method can get satisfactory and fast project results without a large loss of precision.

Navier-Stokes Equations

A mesh with 9,522 triangular real volumes and 4,900 nodes (70-x x 70-h ) was generated. The Reynolds number was estimated based on a flight altitude of 35,000m. The reference speed in this case was taken as the inlet entrance velocity, assuming an entrance Mach number of 15. Hence, at this flight altitude, one can obtain that u_ref = 4,265 m/s. The reference length is 0.1 m, which is consistent with experiments that can be performed in shock tunnels for this configuration. Hence, the Reynolds number can be obtained as

.

The same values of Prd and m were used as in the transonic nozzle and low supersonic airfoil problems.

For this problem, much more severe than the others presented, MacCormack scheme (1985) was not sufficiently robust to guarantee numerical stability. The artificial dissipation terms inherent to MacCormack’s scheme, even with the use of Steger and Warming’s flux vector splitting terms, could not provide the necessary amount of dissipation to initialize the flow simulation. Then, only Jameson and Mavriplis’ scheme (1986) was tested. The CFL number of 0.1 was employed and the following values were adopted for the dissipation coefficients: k⁽²⁾ = 0.65 and k⁽⁴⁾ = 0.04. It is emphasized that the range of variation for values of these coefficients (Jameson, Schmidt and Turkel, 1981) were initially adjusted to transonic flows and not for hypersonic flows. In the last case, the level of required dissipation is much more intense due to the existence of strong shock waves.

Figure 42 shows unstructured mesh used for this viscous case. Figures 43 and 44 show pressure and Mach number contours obtained by Jameson and Mavriplis’ scheme (1986). As a result of this simulation, a boundary layer detachment with formation of a region of circulation was detected. This boundary layer detachment occurred mainly due to interference between upper and lower shock waves. This phenomenon is shown in Figure 45. Figures 46 and 47 show details of the boundary layer detachment at the inlet primary and secondary convergent regions. In terms of convergence history, the results shown for this viscous case were obtained after 50,000 iterations with a reduction of 4 orders of magnitude of the maximum residue. These results indicate that Jameson and Mavriplis’ scheme (1986) presents better robustness features than MacCormack’s scheme (1985) because of the ease with which one can add artificial dissipation for more severe physical problems.

The memory required to storage all main variables and all auxiliary matrices and vectors for both schemes was estimated as: 20 Mbytes for MacCormack’s scheme (1985) and 5.28 Mbytes for Jameson and Mavriplis’ scheme (1986). Such disadvantage of MacCormack’s scheme (1985) is compensated by the high rate of numerical convergence demonstrated in the present work. Table 3 presents a summary of the overall computational results obtained for MacCormack’s (1985) and Jameson and Mavriplis’ (1986) schemes.

Final Remarks

The robustness of MacCormack scheme (1985) for aerospace viscous problems can be improved by implementing a different splitting method as, for example, the van Leer technique (1982). The performance of such method can also be improved using some of the following procedures: the inviscid solution obtained by the Euler equations is stipulated as initial condition; soft start employing a small value for the CFL number and, after some iterations, an increase in its value is processed; and an alternative way of calculating the CFL number, considering flow’s convective and diffusive aspects (MacCormack, 1982, MacCormack, 1985, Mavriplis and Jameson, 1990, and Arnone, Liou and Povinelli, 1991).

Conclusions

Comparisons between MacCormack’s implicit scheme (1985) and Jameson and Mavriplis’ explicit scheme (1986) in the solution of the Euler and Navier-Stokes equations in the two-dimensional space were discussed in this work. Both schemes were tested in three typical problems of the aerospace industry: transonic nozzle, low supersonic airfoil and hypersonic inlet. Their characteristics in relation to computational performance, some aspects of solution quality and robustness were studied, showing advantages and disadvantages between the two algorithms. The robustness of Jameson and Mavriplis’ method (1986) was mainly emphasized in the case of the hypersonic inlet viscous simulation.

In general aspects, the MacCormack’s scheme (1985) is better than the Jameson and Mavriplis’ scheme (1986) in terms of high rate of convergence and run time for the aeronautical problems studied. Its computational cost, associated with memory requirements, is worse than Jameson and Mavriplis’ scheme (1986). General aspects of solution quality and robustness for MacCormack’s scheme (1985) are not so different of those respective features of Jameson and Mavriplis’ scheme (1986). The MacCormack scheme (1985) was not so robust in problems with more severe shock waves. For the hypersonic inlet, the MacCormack scheme (1985) was better than Jameson and Mavriplis’ scheme (1986) in the inviscid case. The meaningful restriction was for the viscous case, where the Steger and Warming’s technique (1981) was not sufficient to guarantee numerical stability. In this case, Jameson and Mavriplis’ scheme (1986) was more robust and able to describe boundary layer detachment. Hence, Jameson and Mavriplis’ scheme (1986) was more efficient than MacCormack’s scheme (1985) for the high speed aerospatial problem studied.

Acknowledgements

This work is part of a Research Project partially financed by Fundação de Amparo à Pesquisa do Estado de São Paulo, FAPESP, through a Doctoral scholarship, level II, to the first author under Process Number 96/2601-4. The authors also thank the partial support of Conselho Nacional de Desenvolvimento Científico e Tecnológico, CNPq, through the Integrated Research Project Grant, under Process Number 522413/96-0.

References

Arnone, A., Liou, M. -S., and Povinelli, L. A., 1991, "Multigrid Calculation of Three-Dimensional Viscous Cascade Flows", AIAA Paper 91-3238-CP.

Azevedo, J. L. F., 1992, "A Finite Difference Method Applied to Internal Axisymmetric Flows", Bulletin of Brazilian Society of Applied and Computational Mathematics, Vol. 3, No. 1, Series II, pp. 1-20.

Batina, J. T., 1990, "Unsteady Euler Airfoil Solutions Using Unstructured Dynamic Meshes", AIAA Journal, Vol. 28, No. 8, pp. 1381-1388.

Chaussee, D. S. and Pulliam, T. H., 1981, "Two-Dimensional Inlet Simulation Using a Diagonal Implicit Algorithm", AIAA Journal, Vol. 19, No. 2, pp. 153-159.

Fox, R. W., and McDonald, A. T., 1988, "Introdução à Mecânica dos Fluidos", Ed. Guanabara Koogan, Rio de Janeiro, RJ, Brazil, 632p.

Harten, A., 1983, "High Resolution Schemes for Hyperbolic Conservation Laws", Journal of Computational Physics, Vol. 49, pp. 357-393.

Harten, A., 1984, "On a Class of High Resolution Total-Variation-Stable Finite-Difference Schemes", SIAMJ. Numerical Analysis, Vol. 21, No. 1, pp. 1-23.

Hefazi, H. and Chen, L. T., 1 992, "A Composite Structured/Unstructured-Mesh Euler Method for Complex Airfoil Shapes", Proceedings of 5th Symposium on Numerical and Physical Aspects of Aerodynamic Flow, Long Beach, California, US, Section 4, pp. 1-6.

Hooker, J. R., Batina, J. T., and Williams, M. H., 1992, "Spatial and Temporal Adaptive Procedures for the Unsteady Aerodynamic Analysis of Airfoils Using Unstructured Meshes", AIAA Paper 92-2694-CP.

Jameson, A., and Mavriplis, D.J., 1986, "Finite Volume Solution of the Two-Dimensional Euler Equations on a Regular Triangular Mesh", AIAA Journal, Vol. 24, No. 4, pp. 611-618.

Jameson, A., Schmidt, W., and Turkel, E., 1981, "Numerical Solution of the Euler Equations by Finite

Volume Methods Using Runge-Kutta Time Stepping Schemes", AIAA Paper 81-1259.

Korzenowski, H., and Azevedo, J. L. F., 1996, "Geração de Malhas Não-Estruturadas para Aplicação em Mecânica dos Fluidos", Proceedings of the 4th Congress of Mechanical Engineering - North/Northeast, Recife, PE, Brazil, pp. 779-783.

Kwon, O. J., and Hah, C., 1995, "Simulation of Three-Dimensional Turbulent Flows on Unstructured Meshes", AIAA Journal, Vol. 33, No. 6, pp. 1081-1089.

Liou, M. -S., 1996, "A Sequel to AUSM: AUSM⁺", Journal of Computational Physics, Vol. 129, pp. 364-382. Long, L. N., Khan, M. N. S., and Sharp, H. T., 1991, "Massively Parallel Three-Dimensional Euler / Navier-Stokes Method", AIAA Journal, Vol. 29, No. 5, pp. 657-666.

MacCormack, R. W., 1982, "A Numerical Method for Solving the Equations of Compressible Viscous Flow", AIAA Journal, Vol. 20, No. 9, pp. 1275-1281.

MacCormack, R. W., 1984, "An Introduction and Review of the Basics of Computational Fluid Dynamics", Lecture Notes, University of Washington, Seattle, Washington, EUA.

MacCormack, R. W., 1985, "Current Status of Numerical Solutions of the Navier-Stokes Equations", AIAA Paper 85-0032.

Maciel, E. S. G., and Azevedo, J. L. F., 1997, "Comparação entre Vários Algoritmos de Fatoração Aproximada na Solução das Equações de Navier-Stokes", Proceedings of the 14th Brazilian Congress of Mechanical Engineering (available in CD-ROM), Bauru, SP, Brazil.

Maciel, E. S. G., and Azevedo, J. L. F., 1998, "Comparação entre Vários Esquemas Implícitos de Fatoração Aproximada na Solução das Equações de Navier-Stokes", RBCM- Journal of the Brazilian Society of Mechanical Sciences, Vol. XX, No. 3, pp. 353-380.

Mavriplis, D. J., and Jameson, A., 1987, "Multigrid Solution of the Euler Equations on Unstructured and Adaptive Meshes", ICASE Report No. 87-53.

Mavriplis, D. J., 1990, "Accurate Multigrid Solution of Euler Equations on Unstructured and Adaptive Meshes", AIAA Journal, Vol. 28, No. 2, pp. 213-221.

Mavriplis, D. J., and Jameson, A, 1990, "Multigrid Solution of the Navier-Stokes Equations on Triangular Meshes", AIAA Journal, Vol. 28, No. 8, pp. 1415-1425.

Mavriplis, D. J., 1995, "An Advancing Front Delaunay Triangulation Algorithm Designed for Robustness", Journal of Computational Physics, Vol. 117, pp. 90-101.

Pirzadeh, S., 1991, "Structured Background Grids for Generation of Unstructured Grids by Advancing Front Method", AIAA Paper 91-3233-CP.

Steger, J. L., and Warming, R. F., 1981, "Flux Vector Splitting of the Inviscid Gas-Dynamic Equations with Application to Finite-Difference Methods", Journal of Computational Physics, Vol. 40, pp. 263-293.

Swanson, R. C., and Radespiel, R., 1991, "Cell Centered and Cell Vertex Multigrid Schemes for the Navier-Stokes Equations", AIAA Journal, Vol. 29, No. 5, pp. 697-703.

Sweby, P.K.,1984, "High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws", SIAMJ. Numerical Analysis, Vol. 21, No. 5, pp. 995-1011.

van Leer, B., 1979, "Towards the Ultimate Conservative Difference Scheme. V. A Second-Order Sequel to Godunov’s Method", Journal of Computational Physics, Vol. 32, pp. 101-136.

van Leer, B., 1982, "Flux-Vector Splitting for the Euler Equations", Proceedings of the 18th International Conference on Numerical Methods in Fluid Dynamics, Lecture Notes in Physics, Vol. 170, Springer-Verlag, Berlin.

Manuscript received: May 1998. Technical Editor: Angela O.Nieckele.

Publication Dates

• Publication in this collection
  29 Oct 2001
• Date of issue
  2001

History

• Received
  May 1998