Numerical Simulation of the Vortex Shedding Behind an Airfoil-Spoiler Configuration

Spoilers are widely used on aircrafts as lateral control devices. Despite their wide usage, very little numerical and theoretical information exist. Numerical simulation of the full unsteady, compressible Euler’s equations over the NACA 23012 airfoil with spoiler is performed using a hybrid Least-Squares finite element/finite difference method coupled to the NewtonRaphson’s linearization technique. The flow patterns behind the spoiler are presented. The pressure coefficient over the upper and lower surfaces are successfully compared to previously published experimental work. The vortex shedding due to the existence of the spoiler is strong specially at high deflection angles. Convection of the vortices will affect the performance of the tail and so a future study of the wing-tail interaction is needed.


Introduction
Spoilers (also called lift dumpers) are devices intended to reduce the lift force of an airfoil in a controlled manner. It was first developed by Martin Aircraft company in (1948). A Spoiler differs from an airbrake in that an airbrake increases drag without affecting lift, while a spoiler reduces lift as well as increasing drag. Spoilers are widely used on aircraft as lateral control devices, as speed breaks, and as lift dampers during landing. Despite their wide usage, very little theoretical information exists. Almost all design work is done by comparison with experimental results followed by wind/water tunnel tests of trail models. The flow field past an airfoil with spoiler is complex. In fact, the flow separated from the upper airfoil surface, due to the adverse pressure gradient generated by the presence of the spoiler, reattaches to the spoiler surface at small angles of attack as predicted by Pfeiffer and Zumwalt [1]. A recirculating bubble is formed upstream of the spoiler hinge. The flow separates again from the spoiler tip and convects into wake as a free shear layer. The flow on the lower airfoil surface also leaves the trailing edge and convects into wake. These shedding vortices make the wake highly turbulent and oscillatory. This unsteady wake affects the effectiveness of the spoiler. Extensive experimental investigations on steady spoiler characteristics have been taken by several authors ( [2]- [5]). On the theoretical side, Abdelrahman [3] developed a numerical scheme based on the steady state incompressible Navier-Stokes equations. Choi et al. [6] studied the transient response of an airfoil to a rapidly deploying spoiler using turbulent compressible Navier-Stokes equations with a turbulence model. Pfeiffer and Zumwalt [1] examined the flow using both water table and oil streaks on a flow-splitter plate mounted on a wind tunnel model. These examinations indicate that the wake has two main characteristics. First it has an unsteady nature due to a shed vortex street. Second, despite this periodic behavior there is a region of wake that remains nearly constant in shape and closes a short distance downstream of the trailing edge. The experimental work of Wentz and Ostowari [2] clearly shows two basic regions: first, an outer essentially potential flow region, and second, a near wake region behind the spoiler. Also, the near wake may be subdivided into two more parts. The part upstream of the wing trailing edge shows very little fluid motion and thus may be considered at constant pressure. The part downstream of the trailing edge is characterized by a pair of vortices; the upper one rotating clockwise (for left to right flow) and smaller one below it is rotating counterclockwise. A numerical study using zonal detached eddy simulation (ZDES) is carried out by Gand [7] shows that the use of zonal detached eddy simulation (ZDES) significantly improves the prediction of the spoiler mean hinge moment compared to the computation performed with a standard Reynolds-averaged Navier-Stokes (RANS)/Spalart Allmaras model on the same grid. Alhawwary et al. [8] simulated numerically the flow field around the BATR airfoil-spoiler configuration with deflected spoiler using the higher order spectral difference method by solving the compressible twodimensional, full Navier--Stokes equations on unstructured quadrilateral grids. Recently Wang and Liu [9] examined the influence of the downward spoiler deflection on the boundary layer flow of a high-lift two-element airfoil consisting of a droop nose, a main wing, a downward deflecting spoiler and a single slotted flap. They found that a downward spoiler deflection results in thicker boundary layer and thicker confluent boundary layer. Also, the flow of spoiler upper surface separates near stall when spoiler deflection is large. The contribution of Geisbauer and Loser [10] give a better insight into the characteristic flow features for applications with deployed control surfaces, for a dynamic spoiler. Exemplary results of an experimental investigation on the steady and unsteady aerodynamic behavior of a configuration with static and dynamic spoiler are presented. Solving the Euler equations is a real challenge due to the convective nature of the equations which necessitates the use of stabilization techniques. Several techniques have been proposed to stabilize the classical finite element method. One of the earliest stabilization techniques is the least-squares (LS) presented by Lynn and Arya [11]. The main advantage of the LS technique is that it produces a symmetric and positive definite coefficient matrix when applied to first order partial differential equations. The main disadvantage of the LS technique is that it produces excessive artificial diffusion in all directions when used with coarse meshes as explained by Taghaddosi et al. [12]. As a remedy for this flaw, the streamline-upwind Petrov/Galerkin (SUPG) method was developed by Brooks and Hughes [13]; in which the diffusion is created along the streamline direction. Then Hughes et al. [14] developed the Galerkin/least-squares (GLS) method which is considered to be a conceptual simplification to the SUPG. The GLS and the SUPG methods have been related to the process of addition and elimination of bubble functions ( [15]- [16]). Brezzi and Russo [17] started the residual-free bubble (RFB) approach, which is further developed by Franca and Russo [18]. The RFB is strongly related to the concept of ''subgrid viscosity'' as illustrated by Guermond [19]. Russo [20] explained why Galerkin method fails when applied to convection-dominated problems. Hughes [21] proposed the variational multiscale method, but Brezzi et al. [22] showed that this method is completely equivalent to the RFB approach. The rest of the paper is organized as follows, in section 2, the balance laws are presented and then linearized. Then in section 3, the finite element formulation is presented followed by the discussion of the boundary condition in section 4. Finally, the findings are presented in section 5.

Governing Equations
For non-polar fluids, conservation of mass, inviscid momentum, and energy, in index notation, with no mass or heat addition, yield the following: The comma denotes partial derivatives. Where ρ is the density,

Linearization and Finite Difference Analysis
One of the most effective linearization techniques is the Newton-Raphson's [12]. In which we set: where H ∆ is the difference between two successive solutions. Discretizing the unsteady term using Euler backward difference and neglecting any higher order terms: where t ∆ is the time step. The linearized form is: where D is a differential operator defined as, (8) f is the right hand side defined as, ( ) n n n n x y I is the identity matrix, and C is given by:

Finite Element Analysis
Since Euler's equations are first order system of partial differential equations, the LS is used since it results in a symmetric and positive definite coefficient matrix which will not be the case if we adopt the SUPG or the GLS techniques. The LS functional to be minimized is defined as: where A is the domain of integration and R is the residual vector, aimed to be zero after the minimization process, given by: Now we can introduce the finite element approximation for the unknowns H  (13) where N is a diagonal matrix containing the isoparametric bi-linear shape functions, e n is the number of nodes per element, and j ∆H is the nodal values. With this approximation; the time dependency is accounted for by the nodal values while the spatial dependency is accounted for by the shape functions. Minimization of the LS functional (11) with respect to the nodal values, H j ∆ , yields local form of the equations as: ∆ = k H r (14) where k is the local influence matrix given by: where a is the element area. Equation (15) shows the element stiffness matrix is both symmetric and positive definite. The right-hand side, r , of equation (14) is given by: Numerical integration, using Gauss-Legendre quadrature, is used to evaluate all the integrations. A home-built FORTRAN code is used to implement the above finite element analysis.

Boundary Conditions
The Euler's equations represent a purely hyperbolic system of equations ( [25]- [26]), and so the number of boundary conditions to be imposed on the domain is determined by the theory of characteristics depending on the incoming/outgoing waves. The details of the boundary conditions that yield a well-posed problem are discussed in [25]. In the finite element method, the numerical boundary conditions are naturally imposed which is considered one of the most important features of the finite element method. For our first order system we have Dirichlet type of boundary conditions except for the slip boundary condition on the airfoil surface. The coordinate rotation method ( [12], [27]- [28]) is used to implement the slip boundary condition. Figure 1 (top left) shows the details of the NACA 23012 airfoil along with the spoiler positioned at 75% of the cord length measured from the leading edge. The thickness of the spoiler is 1% of the cord length and the spoiler length is taken to be 15% of the cord length. Two spoiler deflection angles are considered namely, 30 o and 60 o . The top right section of Figure 1 shows the grid used in the current simulation which is an O-type grid consists of [65x20] bilinear, iso-parametric, quadrilateral elements. The bottom part of Figure 1 is a zoom into the grid around the airfoil as well as a zoom into the grid around the spoiler area. According to the theory of characteristics, a subsonic inlet has three positive eigenvalues and one negative eigenvalue and so one has to impose the characteristic variables correspond to the positive eigenvalues or an equivalent set of the primitive variables. Similarly, for the subsonic exit, one has to impose the characteristic variable corresponds to the negative eigenvalue or an equivalent primitive variable. For the details of the characteristic variables and their primitive variable equivalence, we refer to work of Guaily and Epstein [25]. For an inlet Mach number of 0.2 the following boundary conditions are considered:

Results and Discussion
Upstream: ( , , ) (1, 0.2, 0) u v ρ = Downstream: 0.714 p = Figure 2 shows the evolution of the Mach number contours for a deflection angle of 30 o , using a non-dimensional time step of 0.03 at time stations of 0.6, 1.5, and 2.25, in which the lowest values are shown to be behind the spoiler in the upper part while an increasing behavior is observed for the lower part till the trailing edge where the existence of the spoiler generates wakes as shown in Figure 3 that causes a decrease in the Mach number values. The growth of the wakes behind the spoiler is shown in Figure 3, which is consistent with the experimental observations of Wentz and Ostowari [2]. The evolution of the pressure contours is shown in Figure 4, while Figure 5 shows the distribution of the pressure coefficient over the upper and lower surfaces compared to the experimental results of Abdelrahman [3]. The sudden decrease in the coefficient of pressure to a negative value is due to the presence of the spoiler as described by the contour values in Figure 4. The simulation is repeated for a 60° deflection angle, in which the formation of a recirculating bubble upstream of the spoiler hinge is observed as shown in Figure 6 along with stronger shedding vortices which is consistent with the conclusions reached by Pfeiffer and Zumwalt [1] for large deflection angles.

Conclusion
Numerical simulation of the unsteady flow over the NACA 23012 airfoil with spoiler is performed using a hybrid technique coupled to the Newton-Raphson's linearization method. The usage of Euler-Backward scheme along with Newton-Raphson's linearization technique allows the timestep to be as large as 0.03 which is considerably large when compared to explicit schemes which assures the robustness of the used numerical scheme. The resulting influence matrix from the Least-Squares finite element matrix is symmetric and positive definite which requires only half the computer memory when compared to other techniques like the Galerkin/Least-Squares or Streamline upwind Petrov Galerkin. The iso-parametric bi-linear quadrilateral elements are used in the simulation, but we believe that using subparametric elements e.g. bi-quadratic for the airfoil boundary would enhance the results. The effect of the presence of the spoiler is investigated on the Mach number as well as the pressure fields. Specifically, the coefficient of pressure over the airfoil surface is computed and compared with previously published experimental work with an apparent degree of success. The simulation of the Euler's equations allows for the capturing of the vortices. Specially the one formed upstream of the spoiler hinge at relatively large deflection angles which is successfully captured which is consistent with the literature. The strength of the vortex shedding is studied by considering two deflection angles namely, 30 o and 60 o . The larger the deflection angle of the spoiler the stronger the shedding. This work could be extended to study the effect of the convected vortices from the wing on the performance of the tail. So that the size and location of the tail relative to the wing are chosen in such a way to minimize the negative effects of the conveted vortices.