Sensitivity Analysis of Transonic Flow over J-78 Wings

3D transonic flow over swept and unswept wings with an J-78 airfoil at spanwise sections is studied numerically at negative and vanishing angles of attack. Solutions of the unsteady Reynolds-averaged Navier-Stokes equations are obtained with a finite-volume solver on unstructured meshes. The numerical simulation shows that adverse Mach numbers, at which the lift coefficient is highly sensitive to small perturbations, are larger than those obtained earlier for 2D flow. Due to the larger Mach numbers, there is an onset of self-exciting oscillations of shock waves on the wings. The swept wing exhibits a higher sensitivity to variations of the Mach number than the unswept one.


Introduction
In the 1960s through 1990s, numerical and experimental studies of transonic flow demonstrated formation of double supersonic regions on supercritical airfoils that comprise a flat or nearly flat arc [1,2].Such airfoils may arise in the course of a conventional profile redesign [3] or as solutions of an optimization problem [4].In particular, Jameson [4] designed four airfoils (J-75, J-78, K-80, and GAW-75-06-15) whose thorough studies demonstrated that they admit double supersonic regions and nonunique solutions of the Euler equations in narrow bands of the free-stream Mach number and angle of attack.Richter and Leyland [5] confirmed the nonuniqueness of transonic flow over the J-78 airfoil using an Euler solver with parallelization and adaptation of an unstructured mesh.Hafez and Guo [6] revealed multiple solutions for several symmetric flat-sided and wavy-sided airfoils.Caughey [7] analyzed the stability of inviscid transonic flow over J-78 and Hafez's airfoils using time-dependent boundary conditions.A few more symmetric airfoils supporting nonunique transonic solutions were examined in [8].Meanwhile a nature of the nonuniqueness and the high sensitivity of the flow to variations of boundary conditions were not analyzed.
If a double supersonic region (S-region) forms on an airfoil, then it consists of a bow S-region, which expands as the free-stream Mach number  ∞ increases, and a rear S-region which shrinks and either merges with the bow Sregion or shifts downstream and vanishes.In 2002, Kuzmin [9] noticed that the coalescence of the bow and rear S-regions yields a pressure redistribution on the airfoil and a jump of the lift coefficient.Subsequent numerical studies of inviscid and turbulent flow confirmed the abrupt restructuring of 2D flow field due to the coalescence of S-regions at increasing  ∞ or their rupture at decreasing  ∞ .This phenomenon was scrutinized for a number of symmetric airfoils [10,11] and also for symmetric ones, such as the J-78 airfoil, the Boeing 737 midspan profile [12], and a Whitcomb airfoil with aileron deployments [13].It was shown that both flow regimes (with the double and merged S-regions) can be either steady or exhibit self-exciting oscillations near the trailing edge due to the instability of the separated boundary layer.Figure 1 demonstrates a discontinuous dependence of the lift coefficient on the angle of attack  and  ∞ , for example, for the steady 2D turbulent flow over the J-78 airfoil [11].The lift coefficient is calculated using the standard expression , where  ∞ is the free-stream velocity module,  ∞ is the free-stream density,  is the wing area in planform, and  is the normal force on the wing.
This paper addresses the 3D transonic flow over wings with J-78 airfoil at spanwise sections.The study is focused on the determination of Mach numbers and angles of attack that admit a high sensitivity of the flow field to small perturbations.Such adverse free-stream conditions may occur in International Journal of Aerospace Engineering Figure 1: Lift coefficient   as a function of  ∞ and the angle of attack  for the 2D flow over J-78 airfoil (reproduced from [11] with permission).Sketches on the right illustrate the number of supersonic regions on the airfoil.practice at a vertical gust of wind in cruise flight of civil or transport aircraft.

Formulation of the Problem
We consider the fully turbulent flow past (semi) wings which extend spanwise from the root plane  = 0 to the wing tip plane  = 2. Nondimensional coordinates (), 0 <  < 1 of the J-78 airfoil are given in [11].The trailing edge is normal to the root plane for both unswept and swept wings.In case of the swept wing, the leading edge makes an angle of 10 deg with the normal to the root plane, producing a taper ratio of 0.647346.The airfoil chord length  chord in the root plane is set to 0.5 m.
The wing is inclined at an angle of attack  and enclosed in the box −10 < ,  < 10, 0 <  < 12; see Figure 2. We prescribe a condition of symmetry on the root plane  = 0 and the free-slip condition on the side, top, and bottom boundaries.On the inflow part of the boundary,  = −10, we impose the Mach number  ∞ < 1, the static temperature  ∞ = 250 K, and condition that the velocity vector is parallel to the -axis.On the outflow boundary,  = 10, the static pressure  ∞ = 50,000 N/m 2 is given.The no-slip condition and vanishing heat flux are used on the wing, which is normally assumed to be smooth.In the end of Section 5 we mention also results obtained for a rough wing surface.
The air is treated as a perfect gas whose specific heat at constant pressure is 1004.4J/(kg K) and the ratio of specific heats is 1.4.Then the temperature  ∞ = 250 K determines the sound speed  ∞ = (1004.4× (1.4 − 1) × 250) 1/2 = 316.92m/s.We adopt the value of 28.96 kg/kmol for the molar mass and use the Sutherland formula for the molecular dynamic viscosity.Initial data are parameters of the uniform freestream, in which the turbulence level is set to 1%.

A Numerical Method
Solutions of the unsteady Reynolds-averaged Navier-Stokes equations were obtained with an ANSYS-15 CFX finitevolume solver of the second order accuracy, which is based on a high-resolution scheme for convective terms [15].An implicit backward Euler scheme is employed for the timeaccurate computations.The code utilizes a linearization of the discretized equations and a multigrid accelerated factorization technique for solving the system of linear equations.We use a Shear Stress Transport k- turbulence model, which is known to reasonably predict aerodynamic flows with boundary layer separations [16].
Hybrid computational meshes were constituted by several millions prisms in 35 layers on the wing and over 10 7 tetrahedrons in the remaining region.The nondimensional thickness y + of the first mesh layer on the wing was less than 1.Apart from the boundary layer region, mesh nodes were clustered in vicinities of shock waves and in the wake.
Test computations of transonic flow over the unswept J-78 wing on uniformly refined meshes of approximately 6 × 10 6 , 14 × 10 6 , and 20 × 10 6 cells showed that differences in the lift coefficient obtained on the second and third meshes do not exceed 5%.Global time steps of 2 × 10 −5 s and 4 × 10 −5 s yielded undistinguishable solutions.That is why we employed meshes of 14 × 10 6 cells and the time step of 2 × 10 −5 s for the study of aerodynamic characteristics of the J-78 wings at various  ∞ , .The root-mean-square CFL number (over mesh cells) was about 2.
The solver was verified by computation of transonic flow for a commonly used test case of the ONERA M6 wing at  ∞ = 0.84,  = 3.06 deg.The Reynolds number based on the mean chord is 11.72 million.Figure 3 shows a good agreement of the surface pressure coefficients as computed with the ANSYS-15 CFX solver and a high-order discontinuous Galerkin solver by Bassi et al. [14].
Also verifications of the solver were performed by the simulation of an oscillatory transonic flow past a 18% thick circular-arc airfoil at  = 0,  ∞ = 0.75, and Re = 11 × 10 6 .The obtained amplitude of lift coefficient oscillations was 0.35.This agrees well with the value of 0.37 obtained numerically in [17] using Spalart-Allmaras and Baldwin-Lomax turbulence models.Figure 6 shows maxima and minima of the lift coefficient as functions of the free-stream Mach number.As seen, at Re = 10.4 × 10 6 the amplitude of lift coefficient oscillations is larger than that at Re = 5.2 × 10 6 , though a jump of   is nearly the same.

deg
The leading-edge sweep angle of 10 deg provides a nearly uniform spacing between the bow and rear S-regions on the upper surface of the wing at various spanwise sections.Figure 7 displays the obtained sonic surfaces on the wing at  ∞ = 0.845,  = 0.The height of the bow S-region is smaller than in Figure 4 due to the smaller Mach number, even though the angle of attack is larger than in Section 4.
It can be seen from Figure 7(b) that the height of the Sregions declines with increasing -coordinate.The decline sharpens near the wing tip where 3D effects are essential.
With an increase of  ∞ from 0.845 to 0.849, the bow and rear S-regions approach each other and eventually coalesce in a vicinity of the wing root at  ∞ = 0.849.Contrary to expectations, the coalescence does not produce a jump of the lift coefficient (see the upper solid curves in Figure 8).A rise of   is delayed to the Mach number  ∞ ≈ 0.854 at which the coalescence of S-regions accomplishes over a half of the wing span.Figure 9 shows plots of the pressure coefficient ) on the upper surface of the wing at  ∞ = 0.8535 and  ∞ = 0.8565.We note that the difference of 0.003 in  ∞ means a difference of only 0.95 m/s in the velocity  ∞ .Though the lift coefficient   oscillates in time, the plots of   displayed in Figure 9 are stationary, as the oscillations are associated with static pressure fluctuations on the lower surface of the wing.
When the Mach number increases from 0.857 to 0.861, the lift coefficient declines due to the flow stabilization on the upper surface of the wing and an expansion of the supersonic region beneath the wing.The expansion yields a pressure drop on the lower surface and, hence, the drop of lift coefficient.
At a smaller angle of attack,  = −0.3deg, a dependence of   on  ∞ is similar to the one discussed above for  = 0 though the rise of   is shifted to higher values of  ∞ (see the lower curves in Figure 8).
Dashed curves in Figure 8 demonstrate margins of the lift coefficient oscillations obtained at the vanishing angle of attack for a rough surface of the wing with a sandgrain roughness of 3 × 10 −6 m.As seen, there is a good agreement of the results obtained for the smooth and rough surfaces.

Conclusion
For the swept and unswept wings with J-78 airfoil at spanwise sections and a span of two chord lengths, the numerical study revealed adverse bands of free-stream parameters 0.854 ≤  ∞ ≤ 0.864, −0.6 deg ≤  ≤ 0, in which the lift coefficient exhibits a high sensitivity to small perturbations.In contrast to 2D flow over the airfoil, there are oscillations of the separated boundary layer and shock waves near the trailing edges of the wings.The sensitivity is expected to be even higher for larger wing spans.

Figure 2 :
Figure 2: Sketch of the computational domain.

Figure 3 :
Figure 3: A comparison of the pressure coefficient contours on the upper surface of ONERA M6 wing at  = 3.06 deg,  ∞ = 0.84 as obtained with (a) CFX solver, (b) a high-order discontinuous Galerkin solver [14] (reproduced with permission).
: Lift Coefficient as a Function of  ∞ at  = −0.6 deg When the angle of attack is −0.6 deg, two-dimensional flow over the J-78 airfoil exhibits a lift coefficient jump at  ∞ = 0.84 (see Figure1).As mentioned in Introduction, the jump is caused by a coalescence of S-regions on the upper surface of the airfoil.Meanwhile, computations of the 3D flow over the unswept wing at  ∞ = 0.84,  ∞ = 50,000 N/m 2 , Re = 5.2 × 10 6 showed that S-regions reside at a considerable distance from each other.With the increase of  ∞ from 0.84 to 0.8635 the distance gradually shrinks, while beneath the wing the flow velocity quickly rises.Figure4shows that the S-regions eventually coalesce on the upper surface at  ∞ = 0.864, and this only occurs in the root region.Due to the larger Mach numbers, the terminating shock waves on the upper and lower surfaces of the wing are stronger than those in 2D flow.As a consequence, there is an onset of self-exciting oscillations of the shock waves and separated boundary layer at the rear of the wing.A dependence of the lift coefficient on time, for example, at  ∞ = 0.860, Re = 5.2 × 10 6 , is depicted in Figure5.The oscillations of   are not periodic because of a superposition of different phases of the flow oscillations at different spanwise sections of the wing.

Figure 8 :
Figure 8: Lift coefficient   versus  ∞ for the swept wing at two angles of attack.