A numerical model for thin airfoils in unsteady compressible arbitrary motion

TECHNICAL PAPERS

A numerical model for thin airfoils in unsteady compressible arbitrary motion

Fabiano Hernandes^I; Paulo Afonso de O. Soviero^II

^[e-mail], Empresa Brasileira de Aeronáutica, 12227-901 S. J. Campos, SP. Brazil

^IIEmeritus Member, ABCM, [e-mail], Instituto Tecnológico de Aeronáutica - ITA, 12228-900 S. J. Campos, SP. Brazil

ABSTRACT

A numerical method based on the vortex methodology is presented in order to obtain unsteady solution of the aerodynamic coefficients of a thin airfoil in either compressible subsonic or supersonic flows. The numerical model is created through the profile discretization in uniform segments and the compressible flow vortex singularity is used. The results of the proposed model are presented as the lift and the pressure coefficient along the profile chord as a function of time. The indicial response for the unit step change of angle of attack and unit sharp-edged gust response of the profile are also obtained numerically. The results yielded by the present methodology are also compared with solutions available in the literature.

Keywords: vortex lattice method, compressible flow, unsteady flow; gust; indicial aerodynamics

Introduction

Along the past decade the Generalized Vortex Lattice Method was developed for the unsteady motion, initially for subsonic (Soviero, 1993) and later for supersonic (Soviero and Ribeiro, 1995) and transonic flows (Soviero and Pinto, 2001). In all previous works the profile movement (heaving or pitching) was restricted to the harmonic motion. Therefore, the calculation was performed in the frequency domain.

The only practical way to obtain the aerodynamic loads from arbitrary motions is, according to Bisplinghoff et al. (1955), to use a superposition method (Fourier's integral) of the results obtained for harmonic motions. However, such a methodology is not adequate for sudden movements, which can happen during maneuvers of high performance airplanes, sharp gusts or fast deflections of command surfaces, such as the ailerons. In these cases the number of terms in the series needed to describe forces and moment coefficients can become prohibitively large due to the slow convergence behavior of the solution describing the studied motion.

For the incompressible flow regime, there are well known studies, such as Wagner’s and Küssner’s, that obtained the time evolution of lift of thin profiles with the sudden variation of the angle of attack and profile penetration in sharp-edge gust; both, in fact, indicial responses. In the compressible flow regime, both subsonic and supersonic, a series of indicial responses were presented by Bisplinghoff et al. (1955) as a function of the Mach number for thin profiles. However, calculating these indicial responses analytically is long, tedious, and limited to a short period of time suggesting the need of researching on numeric solutions that are sufficiently fast and applicable to arbitrary motions.

Thus, the study and development of a numeric method that allows the calculation of aerodynamic forces and moments for a profile in arbitrary motion is desirable. Due to the new types of airplanes now in project or construction around the world, in which the aerodynamic efficiency is the main goal, the structures must be lighter and; therefore, more susceptible to, in a disturbed airflow, assuming a motion pattern which may not always be correctly modeled by harmonic or periodic motions. Another application of the indicial formulation is in the helicopter engineering area (Beddoes, 1984).

A numerical model for arbitrary compressible motions, based on the linearized acceleration potential, was developed by Long and Watts (1987). Numerical solutions of the Navier-Stokes equations have been proposed by Shaw and Qin (2000).

A model for two-dimensional incompressible flow for arbitrary motions was developed by Soviero and Lavagna (1997). A similar approach to the one presented here, although limited to the supersonic flow, was developed by Hernandes and Soviero (2002). This model was restricted to the wave equation and presented numerical oscillations in the pressure coefficient results.

An evolution of the work presented by Hernandes and Soviero (2002) is presented here and it considers the classical equation of the non-stationary aerodynamics (the convected wave equation). This new approach solves the problem of the numerical instabilities and extends the range of application of the model for subsonic flows.

Once the indicial response is known, it is possible to obtain solutions for an arbitrary motion using the superposition method employing Duhamel integral, as reported by Bisplinghoff et al., 1955. Lomax et al. (1952) showed how to obtain the flutter speed and stability derivatives also from the indicial response.

The numerical modeling is built through the discretization of a flat plate in uniform segments with compressible vortex singularities. The authors acknowledge that this is the first time this singularity is employed. Lift and pitching moment coefficients, as well as the pressure distribution along the chord as function of time, are obtained from this model. The solutions of the sinking motions and sharp-edge gust response are presented. The solutions are also validated by the solutions available in the literature, such as in Lomax et al. (1952) for the indicial response and Lomax (1954), in Heaslet and Lomax (1947) and in Bisplinghoff et al. (1955) for the sharp-edge gust problem.

Nomenclature

a = sound speed

c_p = pressure coefficient

dt = time interval

n = panel quantity

M = U/a - Mach number

N = time intervals

p = pressure

t = time

t₀ = initial time

U = airflow velocity

U_n = normal velocity over the panel

u' = perturbation velocity in x direction

w = induced velocity in z direction

Greek Symbols

a = profile incidence angle

b = compressibility factor

d f = velocity potential jump

l = sharp-edge gust constant

f = velocity potential

f_x= ¶ f / ¶ x

f _t = ¶ f / ¶ t

f _xx = ¶ ²f / ¶ x²

f_tt = ¶ ²f / ¶ t²

f _xt= ¶ (¶ f / ¶x)/¶ t

r = fluid density

Ñ ² = laplacian operator in x, z domain

Subscripts

j = relative to a panel

k = relative to a time level

Mathematical model

The proposed model is valid for any arbitrary motion, although only two well-known problems available in the literature are addressed here. In the first one, a flat plate with zero incidence angle, immersed in airflow with velocity U, has a sudden change in incidence angle a. Fig. 1 illustrates this indicial response problem. The second one considers the same flat plate with zero incidence angle (immersed in an airflow with velocity U) suddenly exposed to a vertical gust with intensity lU (where l is a constant) traveling against the profile with airflow velocity U. The solution of this problem refers to the response of the profile to a sharp-edge gust entry and is sketched in Fig. 2.

An inviscid fluid is considered. The flow is irrotational and the small disturbances hypothesis is assumed. Therefore, the mathematical model is restricted to the equation of the velocity potential for unsteady flow:

The pressure jump over the chord () is obtained as:

As soon as the angle of attack of the profile is changed to a non-zero angle, a perturbation potential jump df is generated over the profile. In the case of the sharp-edge gust, the potential jumps are generated on the area already reached by the gust. Those potential jumps are immediately substituted by a vortex pair of intensity G and G (where G is numerically identical to df). Fig. 3 illustrates this correspondence:

The intensity of the perturbation potential can be determined by considering the solution for the motion of a piston, which is suddenly placed in an impulsive movement in a compressible flow (Bisplinghoff et al., 1955). The disturbance pressure developed over the profile is given by:

where U_n is the normal velocity over the considered panel. In the time zero, it is equal to Ua for indicial response and to lU (l is a constant) for sharp-edge gust. The variables a and l are assumed unitary in order to facilitate calculations. The results presented are, therefore, function of these quantities. The term of the pressure associated to the impulsive motion is obtained from Eq. 3 and since the term of Eq. 2 related to the impulsive motion is given by the one that includes Df_t , it results in:

Therefore

The velocity U_n , after the initial condition, is given by the initial boundary condition added to the normal velocities induced by the vortices emitted in the previous time levels. The vortices introduced in substitution to the potential jumps are divided in two types: the bound vortices (tied to the profile) and the free vortices (that move with the speed of the undisturbed airflow). The emitted vortices induce normal velocities along the profile chord according to Eqs.5 and 6. The subscript "o" indicates the origin of the vortex (position and time of creation).

where,

X = x - x₀

T = t - t ₀

Numerical model

The proposed model is based on two well-known concepts of theoretical aerodynamics. The first concept is the impulsive generation of vortices in a perfect fluid and the second, the relationship between a pair of counter rotating vortices and the velocity potential jump on the line that links them.

The profile is divided in a convenient number n of panels. The control points, where the velocity potential jumps ( df ) are applied, are placed in the center of each panel and are identified by the index j (l < j < n ). The time levels are identified by the index k (l < k < N ). The variable time step dt corresponds to the elapsed time between iterations.

From the initial motion of the profile, an impulsive motion is generated on each panel, creating a potential jump df, given by Eq.4. Immediately, those n velocity perturbation potential jumps are substituted by pairs of counter rotating vortices.

The sequence of events is explained in Fig.4. Time zero, t = 0 ( k = l ), corresponds to the initial condition, where the piston theory is applied with the boundary condition of normal velocity of intensity U_n over the whole profile chord and posterior substitution of the velocity perturbation potential jumps, df , for pairs of counter rotating vortices. In a next time step ( t = dt , k = 2 ), it is observed that the vortex in the profile trailing edge is free to move with airflow velocity U , (in fact the unsteady Kutta condition) and new velocity perturbation potential jumps, dt , are calculated - now considering, not only the normal velocity U_n , but also the velocities induced by the vortices generated in the previous time t = 0⁺. For further time levels (until k = N ), the assessment of the induced velocities in each panel should consider all the emitted vortices until the previous time level. In the panel edges, a balance of the emitted vortices is made resulting in an algebraic sum in the left extremities of each panel.

The sequence of events for supersonic flow is similar to the subsonic case, just differing in the fact that there are no emissions of free vortices from the profile trailing edge, that is, all the vortices remains tied to the profile. The counter rotating vortices from the perturbation velocity potential jumps for supersonic flow are all bound to the profile because the Kutta condition does not need to be respected in this flow regime. For the subsonic flow case, the vortex of the profile trailing edge (of intensity – and located in the right extremity of the panel j = 1) is, due to model imposition, free to create a wake automatically and to satisfy the Kelvin theorem.

Initially, an explicit integration scheme was studied, where the perturbation velocity potential jumps were calculated as a direct function of the normal velocity over the panel for the considered time level. Moreover, the boundary conditions were established on the profile, i.e., exactly on the control point of each panel. That concept was considered inadequate because the convergence of the method depended directly on the time step adopted. More exactly, when the control point of a considered panel suffered influence from the generated vortex (in the left end of the panel) of that same panel in the previous level, the calculations presented oscillations and divergence in some cases. That divergence was caused by oscillations in the results of the pressure jump coefficient over the profile. The same oscillatory characteristic observed by Hernandes and Soviero (2002) was also noticed by Long and Watts (1987) for the supersonic regime.

Due to these limitations, it was opted not to use an explicit numerical model. A model to establish the boundary conditions halfway between the considered time level and the following one was chosen. Again, that choice was based on the existent analogy between the two-dimensional unsteady flow and three-dimensional steady supersonic flow (Sears, 1954), where the point at which the boundary condition is applied is usually defined in the center of the panel.

Thus, the numerical model becomes implicit because there is, in the considered control point, an influence of its own panel for the considered time level and its solution comes from the solution of a linear system composed by a matrix of coefficients [A] that multiplied by the velocity potential jumps vector, [df]^k, for the considered time level k, results in satisfaction of boundary conditions [W]^k (normal velocities over the panels).

The vector of the perturbation potential jumps is given by:

The matrix [A] is associated with the influence of the generated vortices in a considered time level k and their influence in the same time level. It should still be added to the elements of the main diagonal, the term regarding the impulsive motion, given by 1/(2adt). Thus, the matrix of the influence coefficients is written as:

For subsonic flow, the velocities vector [W]^k is only a function of the normal velocities over the profile due to the profile motion (Ua) added to the velocities induced by the vortices created in the previous time levels. For the method to work out in the supersonic flow, it is essential to consider the additional element of the method - the singularity vortex. In subsonic flows, it is possible to calculate the velocity induced by a vortex in any point of the area affected by its propagation. However, for the supersonic flow, the vortex origin point is singular and it is not possible to calculate the velocity induced in this point. So, it is necessary to define the contribution of the singularity for the velocity field. This concept is well explained in Miranda et al. (1977). The term relative to the induced velocity due to the singularity is defined from the induced velocity in the permanent flow. Eqs.10 and 11 define the resulting velocities vector for subsonic and supersonic flows, respectively.

The system solution is obtained from the equation:

From the solution of the system (vector[dt]^k), it is possible to calculate the aerodynamic coefficients. The circulatory and the non-circulatory (or impulsive) terms of the pressure jump coefficient over the profile are obtained from the equations:

where,

Then, the pressure jump coefficient over the profile results in:

And the lift coefficient for a considered time level is obtained as:

Results and discussion

The aerodynamic coefficients numerically calculated are now presented, for both subsonic and supersonic flows, for unitary change of angle of attack and sharp-edge gust. The calculated results are compared with available solutions in the literature; the results for the indicial response of angle of attack are compared with Lomax et al. (1952) and the results for the sharp-edge gust are compared with Lomax (1954), Heaslet and Lomax (1947), and Bisplinghoff et al. (1955).

The subsonic and supersonic pressure distribution along the chord of a flat plate after an abrupt angle of attack variation are presented, for some time levels, in figures 5 and 6. All results validate the results from literature. In the subsonic case, the pressure distribution tends towards the incompressible and subsonic results, that is, with the characteristic singularity in the leading edge and zero pressure jump at the trailing edge. The supersonic case tends towards the constant pressure distribution of the steady state regime, which is obtained, differently from the subsonic one, after a finite elapsed time.

After pressure coefficient integration along the chord, the lift coefficient for unitary angle of attack is presented in Figures 7 and 8. In the subsonic flow, the steady state values are obtained only asymptotically. The same is not true for the supersonic flow. In fact steady state values are attained earlier as the supersonic Mach number increases. The overall behavior exists due to wave traveling which, in supersonic flow, is only in the downstream direction and in both directions in the subsonic flow. Once again, all calculated values corroborate those from Lomax et al. (1952).

Figures 9 and 10 present the pressure distributions for several time levels for subsonic and supersonic flows unit sharp-edge gusts, respectively. A small oscillation in pressure distribution is noted in Fig. 9 for the two lower time levels and it corresponds to the arrival of the front end of the gust. It is a local behavior due to discretization without influencing the global level of forces.

Time histories of the lift coefficients for subsonic and supersonic unit sharp-edge gusts Mach numbers, respectively, are presented in Figures 11 and 12. Contrary to the problem of step change in angle of attack in both regimes, the lift coefficients start from zero and converge to its steady state value in an infinite or finite time interval depending on the Mach number. All comparisons corroborate the published data.

The last four figures, Figs. 13, 14, 15 and 16, illustrate the variation of center of pressure along the chord with time. In subsonic flow the final position of the center of pressure is at the quarter point from the leading edge. In supersonic flow the final position is, alternatively, at the mid-chord. Both positions in either subsonic or supersonic flow tend toward these steady state values. However, in the initial stages of the motion, there is a remarkable variation in center of pressure position. This fact must be taken into account when simplified versions of the vortex lattice method, where chord discretization is unitary (as in the Weissinger method), are employed.

Finally, it must be pointed out that the results corroborate the references. The largest deviations happen in the abrupt changes of pressure, especially in the supersonic flow, which presents well defined areas in its theoretical curve of the pressure coefficient, as shown in Figs. 6 and 10.

Conclusions

The present work describes a numerical method for arbitrary unsteady motion of thin profiles in linearized compressible flow. The method presents a fast and practical way of obtaining the arbitrary motion response for the subsonic and supersonic flows - an arbitrary motion can be obtained by using a superposition process or simply by changing the boundary conditions since it is a numerical model. Thus, the method proposed is suitable for preliminary aircraft design due to its accuracy, simplicity and computational performance. With this method it is possible to obtain results for a lot of important motions in applied aerodynamics, like gusts and airfoil vortex iteration (AVI) which is very important for the study of the helicopter noise.

It is important to highlight the conclusions concerning the numerical model. It was verified the need to establish the boundary conditions halfway between the considered time level and the following one. This choice is based on the existent analogy between the two-dimensional unsteady flow and three-dimensional steady supersonic flow, where the point at which the boundary condition is applied is usually defined in the center of the panel. This approach, which implies an implicit matrix system, is free from oscillations. Another remarkable feature is the need, in the supersonic flow, to add a term of the induced velocity regarding the singularity vortex.

Initial studies were made with the acoustic approach. That model resulted in oscillations in the pressure jump coefficient (Hernandes and Soviero, 2002). Thus, vortices bound to the profile were chosen, except for the profile trailing edge vortex, which was left free to move downstream in the subsonic case. This model, in spite of being more complex due to the need of circulatory terms calculation, was more appropriate because it did not present the oscillations observed in the acoustic approach mentioned above.

The method yields results that closely follow the available results in the literature (Lomax et. al, 1952, for step change of angle of attack response and Lomax, 1954, Heaslet and Lomax, 1947, and Bisplinghoff et al., 1955, for sharp-edge gust). In particular for supersonic flow, where a significant large number of panels are employed, the error is negligible small, in fact almost imperceptible.

Acknowledgments

The second author gratefully acknowledges the partial support of Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) under the Grant No.302352/2002-3. Special thanks to Edna Mara Marcolino, Luiz Felipe Ribeiro Valentini, Rafael Alves Pereira, and Sheila Regina Braga.

Paper accepted April, 2007.

Technical Editor: Clovis R. Maliska.

• Beddoes, T.S., 1984, "Practical Computation of Unsteady Lift", Vertica, Vol. 8, pp. 55-71.
• Bisplinghoff, R.L., Ashley, H., and Halfman, R.L., 1955, Aeroelasticity, Addison-Wesley, Reading, MA.
• Hernandes, F., and Soviero, P.A.O., 2002,"Numerical Model for Thin Airfoils in Unsteady Supersonic Flow", Proceedings of the 9^th Brazilian Congress of Thermal Engineering and Sciences. Caxambu, Brazil, 7p. (In Portuguese)
• Lomax, H., and Heaslet, M.A., 1947, "Two-Dimensional Unsteady Lift Problems in Supersonic Flight", NACA Report 945.
• Lomax, H., Heaslet, M.A., Fuller, F.B., and Sluder, L., 1952, "Two- and Three-Dimensional Unsteady Lift Problems in High-Speed Flight", NACA Report 1077.
• Lomax, H., 1954, "Lift Developed on Unrestrained Rectangular Wings Entering Gusts at Subsonic and Supersonic Speeds", NACA Report 1162.
• Long, L.N., and Watts, G.A., 1987, "Arbitrary Motion Aerodynamics using an Aeroacoustic Approach", AIAA Journal, Vol. 25, No. 11, pp. 1442-1448.
• Miranda, L.R., Elliot, R.D., Baker, W.M., 1977, "A Generalized Vortex Lattice Method for Subsonic and Supersonic Flow Applications", NASA Contractor Report 2865.
• Sears, W.R., (Ed.), 1954, "High Speed Aerodynamics and Jet Propulsion", Vol. VI, Princeton University Press, Princeton, New Jersey.
• Shaw, T. S., and Qin, N., 2000, "Calculation of Compressible Indicial Response", The Aeronautical Journal, Dec., pp.665-673.
• Soviero, P.A.O., 1993, "Generalized Vortex Lattice Method for Oscillating Thin Airfoil in Subsonic Flow", AIAA Journal, Vol. 31, No. 12, pp. 2380-2382.
• Soviero, P.A.O., and Ribeiro, M.V., 1995,"Panel Method Formulation for Oscillating Airfoils in Supersonic Flow", AIAA Journal, Vol. 33, No. 9, pp. 1659-1666.
• Soviero, P.A.O., and Lavagna, L.G.M., 1997, "A Numerical Model for Airfoils in Unsteady Motion", Journal of the Brazilian Society of Mechanical Sciences, Vol. XIX, No. 3, pp. 332-340.
• Soviero, P.A.O., and Pinto, F.H.L., 2001, "Panel Method Formulation for Oscillating Airfoils in Sonic Flow", Journal of the Brazilian Society of Mechanical Sciences, Vol. 23, No. 4, pp. 401-409.

Publication Dates

• Publication in this collection
  22 Oct 2007
• Date of issue
  Sept 2007

History

• Accepted
  Apr 2007
• Received
  Apr 2007