Unsteady Nonlinear Panel Method with Mixed Boundary Conditions

A new panel method had been developed to account for unsteady nonlinear subsonic flow. Two boundary conditions were used to solve the potential flow about complex configurations of airplanes. Dirichlet boundary condition and Neumann formulation are frequently applied to the configurations that have thick and thin surfaces respectively. Mixed boundary conditions were used in the present work to simulate the connection between thick fuselage and thin wing surfaces. The matrix of linear equations was solved every time step in a marching technique with Kelvin’s theorem for the unsteady wake modeling. To make the method closer to the experimental data, a Nonlinear stripe theory which is based on a two-dimensional viscous-inviscid interaction method for each station along the wing spanwise direction and Prandtle-Glauert rule for compressibility effect were used to enhance the potential results of the method. The fast turnaround time and the ability to model arbitrary geometries is the goal of the present work. Different airplanes configurations were simulated (DLR-F4, light jet, cargo and four engine commercial airplanes). The results of pressure and forces coefficients were compared with the DLR-F4 airplane. The comparisons showed a satisfying agreement with the experimental data. The method is simple and fast as compared with other singularity methods, which may be dependent as a preliminary method to design aircrafts.


INTRODUCTION
Panel methods have been used since mid-1960's to compute the aerodynamic characteristics of complex configurations. The pioneering methods were developed by [1,2]. The methods are presently used as one of the preliminary tools that have been developed to extend its application to the steady, unsteady, viscous flow interaction analyses in subsonic and supersonic flow about complex aircraft, propellers, automobiles, trains, submarines, ship hulls, sails, windmills and other configurations. The panel methods have been used extensively because of their ability to calculate the aerodynamic characteristics in acceptable accuracy at reasonable expense. The discretization is required for the surface of 3D configuration only which is an order of magnitude simpler than volumetric discretization around the configuration needed in finite volume, finite difference, and finite elements methods. Again, for the conceptual and preliminary design studies (less detail is required) this method is superior to give fast and first indication about the aerodynamics analysis of the configuration.
A comprehensive description of the most typical panel methods is presented by [3]. By the 1990's panel methods had largely given way to higher fidelity Navier-Stokes and Euler solvers [4,5]. Despite the promise of Navier-Stokes solvers, accurate viscous drag prediction remained an elusive task. Several researchers started to consider the problem of three-Dimensional Integral Boundary Layer Methods [6,7] with some success. Through 2000's several panel method codes continued the development of higher order approximations to the boundary integral equations [8][9][10][11][12][13][14], however, due to the complexity involved with higher order methods and the lack of robust and efficient integration techniques for higher order approaches, their adoption in the panel method community is limited in comparison with the popular constant collocation type approaches.
Experience with panel methods such as PANAIR, MCAERO, VSAERO, and QUADPAN, developed under NASA contracts and/or industry internal research and development, has shown that low-order methods can provide nearly the same accuracy as higher order methods over a wide range of cases; however, the computation time for low-order panel methods is much shorter than for higher order panel methods [15]. A wide range of lower order and thus more computationally efficient approaches have been developed [16][17][18]. The most recent and most advanced of these approaches is Flightstream, which is based loosely on the surface vorticity solver described in [19]. A variety of applications of panel methods can be found in the literature; see for instance the treatment of unsteady subsonic problems in [20][21][22], transonic and supersonic flows in [23][24] and coupled inviscid boundary layer analyses in [25]. Through the aforementioned previous research, with many of them not mentioned in the current researches related to the panel method and their applications in various fields, it is most likely trying to focus on one method of imposing the boundary condition. Either it was of the Dreichlet method or the Newman condition, or both of them in the case of steady flow only. This means that the users of these methods have limitations to deal with the special shapes that contain zero thicknesses, especially in the Dreichlet method, or those that require more additional calculations due to the nature of the equations therein, as in the Neumann method, and therefore, increase the execution time of the program. Thus, developing a simple method that reduces the execution time and gives some generality to the configurations in calculations of the aerodynamic characteristics is the main motivation for this work and it is a challenge that is not easy. On other side, obtaining the challenge enables solving many calculations that are nonlinear in nature, such as unsteady flow and the accompanying aerodynamic wake behind objects, the effects of the boundary layer, or methods of optimizing the shape, or work of the object and even may be related to fluidstructures interaction to study the aero-elasticity of the plane [26]. Accuracy and efficiency are equal to the complexity that the researcher faces when resorting to use other numerical methods related to as in CFD. The present work, an unsteady flow about airplanes has been considered with mixed surface boundary conditions (Dirichlet and/or Neumann boundary conditions) about three-dimensional complex configurations such as wings fuselage, nacelles etc. in low speed subsonic flow. Prandtle-Glauert rule is used to account for the compressibility effect. Although, the governing equation is linear with Laplace's equation, the problem is still nonlinear because of the position of the vortex sheets wake, which appears as a nonlinearity in boundary condition on the solid surface, as well as in the boundary condi-tion of the vortex sheet itself. As far as forces and moments are concerned, only Navier-Stockes methods are capable to predict the reasonable total drag, so that a viscous effect to the stalling angles of attack is corrected by using a modified nonlinear doublet stripe theory [27]. This method is easy and fast to predict the aerodynamic viscous forces and moments as compared with other methods such as viscous-inviscid method interaction methods, or studying the optimal design of winglets in aircraft wings similar to the research [28], where it will be based completly on the current method with the help of optimization methods to improve the shapes adopted in this field.

MATHEMATICAL FORMULATION
The zero normal velocity boundary condition in an unsteady motion is applied by using body fixed coordinates (x, y, z) as shown in figure 1. where the unsteady motion in other hand is described using stationary coordinate system (X, Y, Z).
By assuming the flow is inviscid irrotational and incompressible for the entire field, then the Laplace's equation can be applied using stationary coordinate system (X, Y, Z); Since the Laplace's equation is not time dependent then the time can be introduced in the boundary condition where the location and orientation can vary with time. Kelvin condition is an additional equation to find the streamwise strengths of vorticity of the wake. The circulation is Γ of the lifting surfaces with its wake are; The kinematic velocity of the surface due to the unsteady motion of the body surface is; where V ∞ is the free stream velocity vector and Ω is the rate of rotation of the body fixed coordinate. Therefore, the instantaneous solution of the Laplace's equation is independent of time derivatives, and steady state solution is used at each time instance, finding the unsteady solution in marching technique. The wake of the lifting surface is the time history of the unsteady motion. Figure 1 also shows a streamlined geometry and its wake. The body surfaces can be represented by S B of the non-lifting bodies and the lifting surface S L with its wake S w . The outer surface S ∞ encloses the flow field of the problem at infinity. Two important regions of interest are contained, the internal and external flow field. The and i are assumed to satisfy the Laplace's equation as in Eq. (1). By applying Green's theorem and combining the external and internal flow regions, the resulting expression yields the following expression for the velocity potential, p , at any point in both regions, which is given by; where r is the distance from point P to the element ds on the surface, and n in the normal unit vector on the body surface S B and S L , which is pointed to the fluid flow. If σ and µ are the source and doublet strength distributions per unit area on the body respectively, then these singularities represent the difference across the boundaries, for the doublet strength it is the potential difference and for source strength it is the normal component difference of the velocity. µ w is the doublet strength distribution over the wake surface S w . The surface S ∞ may be regarded as the far boundary where the conditions only uniform onset flow V ∞ , where the perturbation velocity has essentially disappeared. The velocity potential may be included by ∞ . The Eq. (4) becomes; Alternatively, the velocity of flow field is obtained by gradient Eq. (5) with respect to the body fixed coordinates at the point P. thus; To make the solution unique for the present problem a specification of boundary condition is not enough for the solution, because of two things, the first is the choices of combinations between source and doublet distributions where the source may be used to represent the body thickness effect of the body and doublet is for lifting problems as stated by [3], Secondly, it is the need to fix the circulation around the surface S B which may be expressed by Kutta condition.

Boundary Conditions
Mixed boundary conditions are used to combine direct method (Neumann boundary condition) which specifies zero normal velocity on the solid surface of the body S B . In the present work, the source distribution of the lifting surfaces is equal to zero due to continuous normal velocity through the surface and only doublet distributions (constant doublet lattice method) are used to represent the surfaces of zero thickness. Nacelle of the engine is also treated with doublet panels so that it is only to find the mutual interference between nacelles and other components of the aircraft, where n P is the normal vector at the point P. the source term appears here to consider the possibility of existence of the thick body S B in the geometry.
Indirect method (Dirichlet boundary condition) which is the other boundary condition that specifies a velocity potential everywhere on the boundary surface S B and imposing an internal perturbation velocity potential equals zero for the geometries of non-zero thickness; where σ = -(V ∞ + Ω×r).n P is equivalent to normal.
To make the solution unique, the doublet distributions over the wake surface S w should be known by applying Kutta condition along the trailing edge of the lifting surfaces. Since the wake surface cannot support a load, then the doublet strength distribution is equal to constant along the mean streamlines in wake surface and equals the doublet strength of the trailing edge where the streamline leaves the surface.
where μ T.E. is the doublet strength of the surface panel at the trailing edge. The wake location and doublet distribution are unknown, so that, an iterative method is necessary to solve the problem. More details about the Eqs. (7) and (8) derivation are found in [3]. These expressions cannot be solved directly for the geometry shown in Figure 1, due to the presence of the area integration.

Numerical Solution
The analytical solution of Eqs. (7) and (8) which is stated previously is very hard especially for complex configurations, so that, the problem can be solved numerically with distinct steps as follows, 1. Discretize the geometry using panel method. 2. Forming a matrix of AIC (Aerodynamic Influence Coefficients), which represents linear equations from the Eqs. (7) and (8). The matrix AIC is saved for the next time step. The wake is composed of panels related to the number of time step. 3. Evaluate the solution for these linear equations of the source and doublet distributions over the body surface and wake. 4. Analysis of the results by calculating the velocities and hence the pressure coefficient with forces and moments. 5. Calculate the wake shape using off body velocity equation and repositioned with local velocity flow. 6. Modify the boundary condition at each panel due to influence of new wake shape. For each time step, the upper procedure is repeated, so that, the strengths of source and doublet are changed with time due to the variation of velocity and extending the wake surface with time. The new wake panel travel downstream from the trailing edge in the direction of the local flow velocity with constant doublet strength for each panel at new time step.
The configuration may be divided into two major parts (thick surfaces like fuselage, and/or thin surfaces like wings, tail system and canard, etc. with other components like nacelles and pylon which are treated as thin surfaces.
The trailing edge of the lifting surfaces must be specified for wake surface presentation which Kutta condition satisfied. Each part of the configuration and the wake are discretized into a number of panels of the singularities where previous integrals are now converted into summations. The lifting surfaces are discretized into N L panels and N w panels on the wake. Each panel features doublet distribution whose strengths are constant along the surface of the panel but varying in time. The non-lifting bodies like fuselage, pods, or any other geometry does not generate lifting force which are discretized into the number of panels N B with source /doublet singularities on each panel. The strength of source and doublet are constant along the surface of panel.
The panel has a control point called (collocation point) where Eqs. (7) to (8) (11) and for equation (7) is; The values of the influence coefficients can be obtained in [3]. The panels in Eq. (10) & (12) constitute a set of linear equations with N B + N L doublet strengths as unknowns. All of the wake strengths N W are known except that in the recent region of the trailing edge, which is still defined as a function of the strengths of the doublet panels at the surface trailing edge. These equations must be solved at each time of instance.
The perturbation velocity on the surfaces of S B and S L can be calculated in two ways, the first for the panels on the S B (non-lifting bodies) is calculated by differentiating the doublet strength distribution at any time step as follows; where vx', vy' are the local tangential flow velocities in chordwise and spanwise directions of the panels.
The second is for the panels on S L (lifting surfaces) where the Eq. (6) is used to calculate the local velocities of the flow. By adding the kinematic velocity Eq. (3) to the local calculated velocities then the total local velocity for each panel is used for the time-dependent pressure coefficient as follows; where V T is the total local velocity of flow, ΔCp and Cp are the difference of pressure coefficient for thin surface panel on S L and pressure coefficient on thick panels of S B .
For the low speed subsonic flow, the compressibility effect can be covered by transformation of the coordinate using Prandtl-Glauert rule to reduce the problem to an incompressible flow. This method can be applied at least up to M=0.6 [23].
To include viscous drag and stalling lift phenomenon in the aerodynamic characteristics, results of the wings in particular, the method used in the previous research was modelled to work in the current method [27]. The method is summarized as the calculation of the aerodynamic characteristics by means of using experimental data or a program based on the interaction between the viscous and non-viscous flow for two-dimensional airfoil geometry according to its location along the wing span (stripe location). The angle of attack for that section is adjusted by included the induced angle of the flow and by using iteration with a predetermined accuracy 1x10 5 . One can obtain aerodynamic results that are close to real forces. It must be noted that the method is remodelled using the doublet lattice method instead of vortex ring method. The only difference is the position of the vortex rings and control points.   [29]. The processing of data for the body and the wing was published in [30]. As shown in the figure, the body consists of sections distributed along its length from the nose to the tail. These sections consist of points distributed circumferentially in each section. The sinusoidal distribution method of these sections was used to condense them in an area where the curvature is near the nose and tail of the body. The aircraft models have been presented in Sumo 2.7.10 program which is a 3-D graphical tool for various complex aircraft shapes. Also, it is used for preparing the domain of consideration to be used in computational fluid dynamics. The wing and body were completely modelled by using this program. The total number of panels which represents the aircraft is 3500 panels. The chord wise divisions between the wing and fuselage have been taken with same divisions at the connection region which makes the region to be easy to treat without having to divide the upper part and the lower part of the wing and body as shown in Figure 3. Its location, orientation, and dihedral angle are also fixed separately (not included in the airfoil data) for easy translation to any position along the body. The connection line between the fuselage and the wing is used to represent the entire fuselage panels. The fuselage in this case is divided into three parts, as shown in the figure. The central part, which represents the wing's connection to the fuselage whose sections are re-divided according to the points distributed along the line of contact.
This, in turn, requires that the front and back parts of the body of plane have the same subdivisions of the circumferential points of the center part of the plane. It is possible to make more than one place of contact between the wing and the fuselage, for example, tail, wing and canard in the same way discussed before. Since the shape is symmetrical around the x-z plane the other side has been copied to perform calculations on the entire plane model. The distribution of panels for the body and wing was considered by the trial and error. The results indicate that the optimum division is achieved by having a mean panel aspect ratio of 0.3. To make the plane have more than one component, different models of airplanes were simulated to check the program's ability to deal with complex shapes.
All cross-sections airfoils of the wing have been taken to be a half-thickness line. This will reduce a lot of panel's number and effort during wing panelling. Figure 4. shows planes with winglets and nacelles. As evident from the figure, the possibility of building complex forms of aircraft with different bodies and shapes is an easy job.
From aforementioned aircraft models, the possibilities of building aircraft using this method are greater and faster than the methods used to construct models with thick wings. There is no need for special treatment for the upper and lower surfaces of wing to the connection region between wings and bodies. The program will be developed to treat the remaining parts of the plane, such as the pylon, wheels, etc.

RESULTS AND DISCUSSIONS
In this section, the most important results of the present method from the aerodynamic point of view are discussed and compared with the experimental results of DLR-F4 airplane obtained from [28 and 29], which give the present work reliability in analysing aircraft. This section also illustrates the applicability of the current method to find the mutual interference between the fuselage, wing, nacelles, etc. In the present results, here are two types in presentation of pressure coefficients distribution on the aircraft surfaces. The first are the distribution of the pressure coefficient C P on the parts that contain thicknesses such as the fuselage or stores under the wing, etc. These bodies are treated using the Dirichlet's boundary condition, which calculates the value of pressure coefficient on the surface directly. The second presentation displays the value of the difference in the pressure coefficient ∆C P for the lifting surfaces which are treated by using the Neumann's boundary condition. Therefore, the minimum and maximum values of the pressure coefficients for each type differ according to the type of the boundary condition on the surfaces.
The red colour represents the maximum value of pressure coefficient, passing through yellow and green colours, while the blue colour represents the minimum value. Figure 5 shows the distribution of the pressure coefficient on the model of DLR-F4 wing-body configuration. It is clear that the pressure coefficient especially at the contact line between wing and fuselage of the aircraft shows a smooth distribution. The continuous distribution of doublet panels between the body and the wing prevent a strong vortex at the line of contact, although the boundary conditions of the two surfaces are different. As shown in the figure, the effect of the wing on the fuselage pressure coefficient is evident on the panels surrounding the contact line of contact the wing and the fuselage.
The maximum values for the distribution of the pressure coefficient of the fuselage are the stagnation areas, which are clearly shown near the nose and tail of the airplane. On other hand, the wing shows a red colour of ∆C P at leading edge of the wing, where the differences between the upper and lower surfaces are high. For the trailing edge the difference is nearly equal to zero. The current method appears to be successful in predicting various aerodynamic interaction effects where the regions of high and low pressure are clearly illustrated in the figure.
For the sectional pressure, Figure 6 shows the wing sectional ∆C P of the wing-body configuration. Although the effect of wing thickness has been cancelled, it seems that the distribution of ∆C P along the wing chord-wise direction is in fair agreement with experimental data. There is a slight increase in the pressure difference values at the leading edge of the wing; this is because of the simple theory used in the work which depends only on the potential theory, but not on the treatment to the suctional vortex at the leading edge of the wing. It must be mentioned that, the viscous effect is included only when evaluating the total airplane forces and moments. All experimental data on the span wise positions of the wing weretaken to be the ∆C P for the lower and upper surfaces pressure coefficient of the wing. Figure 7 shows the relationship between the lift coefficient and the angles of attack of the aircraft. The agreement between theoretical results and the experimental data is very good as shown because of the effect of flow viscosity with a non-linear stripe's theory, which enforces aerodynamic characteristics as that of two-dimensional viscous flow data where the separation had been included. Although the compressibility effect was calculated using Prandtl-Glauert correction method, the results were close to the experimental data.  Figure 8 shows the relationship between drag coefficient and angle of attack of the aircraft. The program seems to be well estimated as compared with the experimental data. Regardless that there is a discrepancy between the theoretical and experimental results at high angles of attack of the aircraft due to no consideration of the viscous drag of the fuselage but, generally, it gave very good results to solve a complex model such as the current one. For the pitching moment coefficient, Figure 10 shows a good agreement especially for the positive angles of attack.
Some of the points are deviated from the experimental data, which may be attributed to high separating boundary layer on the lifting and non-lifting surfaces. The wake shape for the current aircraft model as shown in Figure 11, where it clearly shows the growth of the trailing vortices at wing tips, also the starting vortex resulting from the sudden acceleration of the aircraft. The wake surface shape near the fuselage region has small twisting shape as that in the wing tips, which indicates that the method successfully treat the connection between the wing and the fuselage where no strong vortex is created, which was caused by a discontinuous of bound vortex along the span of the wing. Finally, the most important note in the wake field is the size of panels, which are similar to the element size on the wing, and this is necessary to give a good result of accuracy. A simple equation can be generated to give the time value used to determine the length of the panel along the wake as shown in Figure 2. this equation is; where Δt is the time step in the unsteady motion calculations; b is the aircraft wing span, F b represents a factor multiply by the wing span which gives the total length of the wake surface usually for good precession greater than 2.0. N w is the number of wake panels along each stripe. To study different models of aircraft in terms of its wing position and shape as well as its wing connection to the body of the plane. Figure 12 shows pressure coefficient distribution of the models of light jet, cargo and four engine transport aircrafts at angle of attack 3 o and M ∞ =0.16, Re=3x10 6 . The number of panels for these aircraft are (2420, 2968 and 3664 panels) respectively. The presence of the winglet, nacelles gave the program a great flexibility to achieve complex aircraft configurations, of whatever sizes, small or large. The mutual interference effects between the engine nacelle and the fuselage, wing (especially at the location where are close to each other) are clearly illustrated in these figures. Thrust force is not included in calculation, the only aerodynamic interference effects of the engine's casing (nacelle) are studied instead. Figure 13 shows the forces lift and drag coefficients as a function of the angle of attack. The results show that the four-engine aircraft had a greater drag coefficient as compared with other types and stalling earlier for Mach no. (0.16). The wake consideration for these aircraft is shown in Figure  14. The wakes representations behind these models of aircraft show the principle of interaction between different parts of the airplane which formed this wake. The presence of the wake of the wing and engines led to intense interaction between them. On the other hand, the effect of calculating the boundary layer is to reduce the rate of error in the results, as noted in the DLR-F4 results when compared with experimental data. The efficiency of the present method for the computational running time is more than of other methods which calculate the domain of consideration like FVM of FDM or even the singularities methods that solve the fully thick bodies. Also, the results show that the accuracy is satisfactory as compared with experimental data. The computing time for DLR-F4 airplane to find the aerodynamic characteristics is approximately 2 min in Core i7 sixth generation desktop personal computer. The program is built completely using MATLAB 2012a.

CONCLUSIONS
A nonlinear unsteady panel method program had been built using the Matlab 2012a environment to calculate the aerodynamic characteristics of different types of aircraft, like DLR-F4, cargo, light jet, four engine commercial airplanes. The mathematical work was based on the first-order singularity method with mixed types of boundary conditions (Newman for lifting surfaces and the Dirichlet for thick bodies). The Prandtl-Glauert rule was used to compute the compressibility effect of the flow, while stripe method was used to calculate the stalling criterion and drag force of the airplane. Theoretical results obtained from the present program were shown and compared with published experimental work of DLR-F4 airplane. The following conclusions were obtained: 1. The current method represents an easy method to deal with steady and unsteady flow especially for the connection area between lifting and other parts of airplane.
2. This method takes number of panels less than that in other methods, especially those with thickening in the lifting surfaces.
3. The program can calculate a nonlinear force for the steady or unsteady flow with the interaction between the parts of the airplane and the wake. 4. The experimental results were compared at different attack angles to obtain reliability of the present results for the DLR-F4 airplane, and good agreement between them was shown. 5. A comparison was made with the rest of the other types of aircraft and showed that it gives quick results which may be used to optimize the airplane configuration as a future work. Also, the existence of the stripe theory succeeded the predicting results of the stalling angle and drag force. It is possible to compute many examples related to wind turbines, exchange interactions between different parts, and wake rollup behind the airplane. The program is currently developed for the modelling of pylon and thrust force of the engines.