504 IIUM Engineering Journal, Vol. 15, No. 1, 2014 Sulaeman and Ahmed 71 APPROXIMATE FUNCTION FOR UNSTEADY AERODYNAMIC KERNEL FUNCTION OF AEROELASTIC LIFTING SURFACES ERWIN SULAEMAN AND LAYEEQ AHMED Department of Mechanical Engineering, Faculty of Engineering, International Islamic University Malaysia, PO Box 10, 50728, Kuala Lumpur, Malaysia. esulaeman@iium.edu.my ABSTRACT: Prediction of unsteady aerodynamic loads is still the most challenging tasks in flutter aeroelastic analysis. Generally, the numerical estimation of steady and unsteady aerodynamics of thin lifting surface is conducted based on an integral equation relating aerodynamic pressure and normal wash velocity. The present work attempts to increase the accuracy of the prediction by using an approximate approach to evaluate kernel function occurring in the integral equation in the form of cylindrical function. Following previous approximation approaches by other researchers to solve the cylindrical function for planar lifting surfaces, this paper extends such approaches to non planar lifting surfaces. To increase the accuracy of the method, the integration region of the kernel function is divided into two parts – near and far regions, where a nonlinear regression curve fitting technique is adapted to estimate the denominator part of the cylindrical function of each region. ABSTRAK: Penelahan daya aerodinamik tidak stabil merupakan satu tugas yang mencabar dalam menganalisis getaran aeroanjalan. Umumnya, anggaran berangka untuk daya aerodinamik stabil dan tidak stabil pada permukaan mengangkat yang nipis, adalah berdasarkan kepada persamaan kamiran di antara tekanan aerodinamik dan halaju aliran udara pada garis normal yang terhasil di bawah sayap pesawat. Kajian ini adalah bertujuan untuk menghasilkan penelahan daya aerodinamik yang lebih tepat dengan menggunakan pendekatan kira hampir untuk menilai fungsi Kernel yang terdapat dalam persamaan kamiran dalam bentuk fungsi silinder. Dengan menggunakan pendekatan kira hampir yang digunakan oleh penyelidik sebelumnya untuk menyelesaikan fungsi silinder pada permukaan mengangkat satah, kajian ini mengembangkan pendekatan tersebut kepada permukaan mengangkat tak sesatah. Untuk meningkatkan lagi ketepatan penelahan, kawasan pengamiran fungsi Kernel dibahagikan kepada dua bahagian, kawasan hampir dan kawasan jauh, di mana penyuaian lengkung regresi tak linear digunakan untuk kiraan hampir penyebut pada fungsi silinder pada setiap kawasan. KEYWORDS: aeroelasticity; unsteady aerodynamics; kernel function; cylindrical function; lifting surface 1. INTRODUCTION Investigations are in progress for the improvement of accuracy and efficiency of the load prediction methods for countless types of aerodynamic configurations ranging from simple two dimensional airfoils to complicated full scale aircrafts [1]. Since its derivation in 1940, Küssner’s governing formulation is the origin of most aerodynamic load formulations for thin lifting surfaces [2,3]. The fundamental part in the Küssner integral equation is the so-called kernel function which relates aerodynamic pressure and normal wash velocity. Several forms of the kernel function have been proposed in the past, including the formulation of Watkins et al. [4], Laschka [5], Yates [6], Landahl [7], and IIUM Engineering Journal, Vol. 15, No. 1, 2014 Berman et al. [8]. The effect of acoustic [9]. The widely accepted formulation of the kernel function is due to Landahl where the formula is used in aeroelastic analysis tool of MSC formulations contain a hyper whose solutions are not readily available in commercial software. Most of the literatures use approximation methods to solve the incomplet function [3, 8]. The first approximate method was adopted by Watkins et al. used in unsteady aerodynamic doublet lattice method of Ref. 8. Laschka both analytical and numerical Similar numerical approximation of series of fraction with constant numerators was offered by Dat and Malfois [13,14] in order to minimize the error in the approximation. Ueda expansion series to solve analytically the problem. However, the Ueda series requires a large number of steps for a large number of arguments of oscillating fun [2]. To increase the accuracy of the kernel function evaluation, an analytical separation of singular and regular functions occurring in the incomplete cylindrical function is proposed in [16] by modifying the expansion series. Another a Bismarck-Nasr [17] based on extension of the kernel function formulation to transonic flow and separated flow. and Bliss [2] suggested an approximation to th dividing the regions of integration into two parts which are near and far fields. Their work in Ref. 2 is presented for the planar lifting surfaces. Following the Epstein in the present work such approac accuracy is increased by adding more terms to the approximation functio Fig. 2. GOVERNING EQUATION The Küssner integral equation relating in the unsteady potential derivation is written as [3]: ol. 15, No. 1, 2014 Sulaeman 72 l. [8]. The effect of acoustics in the integral equation is formulated by Yu et al. The widely accepted formulation of the kernel function is due to Landahl where the eroelastic analysis tool of MSC Nastran software. All of these formulations contain a hyper-geometric function of incomplete cylindrical function type solutions are not readily available in commercial software. Most of the literatures use approximation methods to solve the incomplet The first approximate method was adopted by Watkins et al. used in unsteady aerodynamic doublet lattice method of Ref. 8. Laschka [11 both analytical and numerical solutions for the kernel function with three digit accuracy. Similar numerical approximation of series of fraction with constant numerators was offered by Dat and Malfois [12]. Least square technique is suggested by Desmarais in order to minimize the error in the approximation. Ueda [15] expansion series to solve analytically the problem. However, the Ueda series requires a large number of steps for a large number of arguments of oscillating functions as shown in [2]. To increase the accuracy of the kernel function evaluation, an analytical separation of singular and regular functions occurring in the incomplete cylindrical function is proposed by modifying the expansion series. Another analytical solution is presented by based on differential equation approach. [18] described possible extension of the kernel function formulation to transonic flow and separated flow. and Bliss [2] suggested an approximation to the incomplete cylindrical function by dividing the regions of integration into two parts which are near and far fields. Their work in Ref. 2 is presented for the planar lifting surfaces. Following the Epstein-Bliss approach, in the present work such approach is extended to non planar lifting surfaces and the accuracy is increased by adding more terms to the approximation function. Fig. 1: Surface and Panel (box) geometry. GOVERNING EQUATION The Küssner integral equation relating the pressure and the normal wash distribution in the unsteady potential derivation is written as [3]: Sulaeman and Ahmed the integral equation is formulated by Yu et al. The widely accepted formulation of the kernel function is due to Landahl where the Nastran software. All of these geometric function of incomplete cylindrical function type Most of the literatures use approximation methods to solve the incomplete cylindrical The first approximate method was adopted by Watkins et al. [10] which is 11] presented solutions for the kernel function with three digit accuracy. Similar numerical approximation of series of fraction with constant numerators was . Least square technique is suggested by Desmarais presented an expansion series to solve analytically the problem. However, the Ueda series requires a ctions as shown in [2]. To increase the accuracy of the kernel function evaluation, an analytical separation of singular and regular functions occurring in the incomplete cylindrical function is proposed nalytical solution is presented by described possible extension of the kernel function formulation to transonic flow and separated flow. Epstein e incomplete cylindrical function by dividing the regions of integration into two parts which are near and far fields. Their work Bliss approach, h is extended to non planar lifting surfaces and the the pressure and the normal wash distribution IIUM Engineering Journal, Vol. 15, No. 1, 2014 Sulaeman and Ahmed 73 �� � � � ∆� �, , �, �, �, �, �, �� 8� �� � � �� �� (1) where w is the normal wash velocity at the control point (x,y,z) on the lifting surface as shown in Fig. 1, U is the free stream velocity of the unperturbated flow assumed in the x axis direction; and �� is the nonstationary aerodynamic pressure difference at point � , � , �� along the line of doublets of each trapezoidal aerodynamic model box shown in Fig. 1. The free stream dynamic pressure is denoted by q and is equal to ρU 2 /2 where ρ is the free stream air density. The kernel function of the integral K is a function of relative position, the free stream Mach number M, and the reduced frequency k. The relative position between the control point and the doublet pressure location can be expressed as: �� � � � � � � � � �� � � � � and the reduced frequency k is defined as � � � � where L is the reference length. 2.1 Planar Surfaces The kernel function K given in Epstein[2] for planar lifting surface is written as ! �� , �, �, �" � # $%&'( �# $%&) *√,� - � � - . �, �, ,� (2) where B is the incomplete cylindrical function formulated according to Landahl [7] and Ueda [15] as . �, �, �� � / # %&0 � � - 1��2� ) $3 �1 � 1� � / #%&56 1 - 7��2� �7 ∞ $)5 (3) and the modified distance r , R and X are defined as � � 8 � ��� � | �| * � :��� - ;�� � , � �� � �*;� � � ;� � 1 � �� Separating the real and imaginary terms [2] from Eqn (3) yields . �, �, ,� � .5 - < .% � 1�� / cos ��7� 1 - 7��2� ∞ $)5 - <� � / sin ��7� 1 - 7��2� ∞ $)5 (4) IIUM Engineering Journal, Vol. 15, No. 1, 2014 Sulaeman and Ahmed 74 and integrating the equation by parts gives . �, �, ,� � 1� � !V"$|)/5|∞ � ��� � / # %&56 D 7��7 ∞ $)5 (5) where the term V represents the non-integral term E � #%&56 � 7√1 - 7� � 1� � ��# $%&56 F�7 - 81 - 7�G (6) and the function F(u) is defined by D 7� � �7 - 81 - 7� (7) In the present work, a simple curve fitting technique is utilized to approximate F(u) as follows: D 7� H I 7� � J IK 7� for 7 N 1.5 IQ 7� for 7 R 1.5 S (8) where fa (u) and fb (u) are defined as IK 7� � T U% 7% V %W� (9a) IQ 7� � T X% 7$% Y %WZ (9b) and the coefficients ai and bi are defined in Table 1. Table 1: Coefficients ai and bi of Eqn (8) and (9). i ai bi 0 0.99999933 -- 1 -0.99968242 0.49998705 2 0.49488907 0.00017800315 3 0.0303999182 -0.12594842 4 -0.21413252 0.00099707005 5 0.13894299 0.071073596 6 -0.041028643 -0.032656364 7 0.0048280266 -- The accuracy of the present approximation is demonstrated in Fig. 2 and Fig. 3. In Fig. 2, the plot of f(u) resembles the target function F(u). The plot of the error representing the difference between F(u) and f(u) is presented in Fig. 3. The maximum error for the present approximation is 3.5 x 10 -6 at u = 1.4. Figure 3 also presents the plot of the error if the approximation of Epstein-Bliss [2] is used. The maximum error of their error is 0.015 which occurs at u = 0.5. Therefore, the present approximation provides much better accuracy compare to the approximation of Ref. 2. Note that the plot for the Epstein-Bliss IIUM Engineering Journal, Vol. 15, No. 1, 2014 Sulaeman and Ahmed 75 approximation in Fig. 3 is multiplied by 10 -3 whereas for the present approach is multiplied by 10 -6 in order to show the sensitivity changes in the error. Fig. 2: Comparison between the present approach f(u) and the target function F(u). Fig. 3: Comparison between the error of the present approach f(u) and the Epstein-Biss approach with respect to the target function. Figure 4 shows comparison of Br(k,r,X) between the exact analytical solution of series expansion [14, 15] and the present approximation. Figure 5 shows comparison of Bi(k,r,X) between the exact and the present approximation. To show the accuracy of the present approximation, the error or the difference between the analytical solution and the present approximation is plotted in Figs. 6 and 7 for Br and Bi respectively. Both figures also present the error for Br and Bi calculated using Epstein-Bliss approximation which demonstrates the improvement of the present approximation. Figure 6 shows that the maximum error of the Epstein-Bliss approximation for Br is 0.0127 whereas the present approximation gives the maximum error of 0.00047. In Fig. 7, the maximum difference for Bi using the Epstein-Bliss approximation is 0.0094. But the new approximation plot shows the accuracy of 0.0000667 which increases the applicability of the new approximation to a good extent in solving the incomplete cylindrical function in the kernel function. 0 0.2 0.4 0.6 0.8 1 1.2 0 1 2 3 4 5 F (u ) u Exact Present Approx -20 -15 -10 -5 0 5 0 2 4 6 8 10 E rr o r (E x a c t - A p p ro x ) u Epstein-Bliss ( x 10-3) Present Approx. ( x 10-6) IIUM Engineering Journal, Vol. 15, No. 1, 2014 Sulaeman and Ahmed 76 2.2 Nonplanar Surfaces Rodemich [19] derived an expression in 1965 for the kernel function of nonplanar lifting surfaces as follow: � # $[%\'(] ^ Z_Z - �_��/�Z� (10) where _Z � `ab cd � c5 � _� � 1�Z� �� `ab cd � � b