Flutter Analysis of a 3-D Box-Wing Aircraft Configuration

The aim of the current paper is to present a flutter analysis of a 3-D Box-Wing Aircraft (BWA) configuration. The box wing structure is considered as consisting of two wings (front and rear wings) connected with a winglet. Plunge and pitch motions are considered for each wing and the winglet is modeled by a longitudinal spring. In order to exert the effect of the wing-joint interactions (bending and torsion coupling), two ends of the spring are located on the gravity centers of the wings tip sections. Wag ner unsteady model is used to simulate the aerodynamic force and moment on the wing. The governing equations are extracted via Hamilton’s variational principle. To transform the resulting partial integro-differential governing equat ions into a set of ordinary differential equations, the assumed modes method is utilized. In order to confirm the aerodynamic model, the flutter results of a clean wing are compared and validated with the previously published results. Also, for the validation, the 3-D box wing aircraft configuration flutter results are compared with MSC NASTRAN software and good agreement is observed. The effects of design parameters such as the winglet tension stiffness, the wing sweep and dihedral angles, and the aircraft altitude on the flutter velocity and frequency are investigated. The results reveal


Introduction
Box Wing Aircraft (BWA) is one of innovative configurations commonly used by airplane designers to reduce the aircraft emissions. Initial investigations on the box wings were implemented by Prandtl 1 who introduced the box wing aircraft as the best wing system that dramatically decreases the induced drag. A reduction of induced drag in the BWA has significant influence on the airplane weight and performance. Moreover, the reduction of induced drag has other benefits such as reduction of geometric dimensions and fuel consumption. On the other hand, the reduction of the wing span, due to using front and rear wings, could decrease the traffic problems at airports for midsize and large airplanes such as Boeing 747, Airbus 380 and Antonov 225.
Due to the great importance of the box wings, many researches have been conducted in the field of design and consequently, various types of configurations have been introduced to the international scientific society such as diamond wings, struts-and truss-braced wings, box-wing and PrandtlPlane configurations. Wolkovitch, 2 Livne and Weisshaar, 3,4 and Cavallaro and Demasi 5 comprehensively reviewed the previous efforts on the joinedwing configurations, and accurately expressed the challenges, ideas and innovations on the configurations. In a Master of Science thesis at Pisa University, the basic design of a 250 passenger PrandtlPlane aircraft (PrP250) was investigated. 6 Torenbeek et al. 7 have discussed the important challenges of European Union such as the quality of aircraft wings, environmental issues, flight safety and air fleet performance in this work. In addition, he compared the induced drag between an optimal box wing and conventional wing. Frediani et al. 8 presented the methodologies used to design a class of PrandtlPlane ULM aircraft. The methodologies met the design requirements such as high visibility, low noise in cabin, and low stall speeds and enhanced the safety. Jansen and Perez 9 inquired the multidisciplinary design and optimization of non-planar configurations. They considered the coupling between aerodynamic and structure for different aircraft sizes and mission requirements. The results revealed the performance effects on the design of the configurations. Barcala et al. 10 stated that the existence of the front and rear wings in different planes can make the aerodynamic analysis of box wing configurations by theoretical/computational methods even more difficult. Therefore, they analyzed an experimental investigation on several configurations of the box wings. In the research project "IDINTOS" at Pisa University, Frediani et al. 11 designed and manufactured an ultralight amphibious PrandtlPlane prototype and some of the benefit results were figured out by means of both numerical and experimental analyses. There are also some other efforts in the field of joined-wing design which have evaluated the design requirements of the different configurations such as aerodynamics, flight mechanics, and weight optimization. [12][13][14][15][16] Although numerous studies have been conducted on wing aeroelasticity, 17-25 a preliminary literature review shows that most previous studies on box wings are focused on design and aerodynamic performance. Lange et al. 26 studied the design of a transonic box wing and remarked that "because of the complexity of the configuration, proper assessment of airplane characteristics utilizing this concept can only be accomplished by the use of rather elaborate analytical procedures for aerodynamic, structural loading and aeroelastic characteristics". Moreover, the existence of two wings separated in different planes has made the BWA aeroelastic analysis by theoretical or computational methods a difficult task. Hence, the vast majority of researchers have recently investigated the BWA aeroelastic analysis via professional software, [27][28][29][30][31][32] and a limited progress has been made on the mathematical modeling of aeroelastic BWA configurations.
Regarding the nonlinear aeroelasticity of the joined-wing configurations, in the recent years a number of studies have been published employing a discrete time-domain statespace approach which can convert the nonlinear system to a set of sub-linear systems. Each sub-linear system is considered to be linear whose response can be computed by a straight forward time integration. These studies are focused on the nonlinear systems due to large deformations (Chao et al. 33 ), geometric nonlinearities (Blair et al., 34 ) and nonlinear aerodynamic effects (Cavallaro et al. 35 and Alizadeh et al. 36 ).
In the linear aeroelastic analysis, flutter speed and frequency are obtained with a smaller value compared to nonlinear cases. If the wing aspect ratio is low, the linear analysis can be accurate. On the other hand, because of the simplicity of solving the linear governing equations, many researchers prefer to utilize the linear aeroelastic analysis. Mirabbashi et al. 37 analyzed the flutter of an airfoil carrying a flexibly mounted unbalanced engine by analytical and experimental approaches. The typical wing section and the connected engine are modeled by a 5DOF system. The results revealed that some parameters of the engine can decrease the flutter boundary e.g. mass, pylon length and thrust. Ajaj et al. 38 studied the flutter of a telescopic, multi-segment and span morphing wing. They considered the wing as a stepped and cantilever Euler-Bernoulli beam and showed that the span morphing can augment the flutter boundary.
Durham and Ricketts 39 concluded that the particular joined-wing configuration can have higher flutter velocity in the high altitude compared to high aspect ratio configurations and they could prove their claim using wind tunnel tests. Van Aken 40 perused the feasibility of using active controls to postpone the whirl-flutter on a joined-wing configuration. He utilized the CAMRAD/JA code to extract the governing equations (a set of linear differential equations). The flutter analysis showed that the active control can delay the onset of tilt-rotor whirl-flutter. Lee 41 developed a methodology for aeroelastic design of the joined-wing configuration. The stability boundaries were obtained using Rayleigh-Ritz and finite element methods and the effective parameters on the aeroelastic behavior were also assessed. The results illustrated that the body-freedom flutter velocities decrease as fuselage flexibility increases. Furthermore, the center of gravity location and pitch moment of inertia can affect flutter boundary. Canto et al. 42 expressed the PrandtlPlane aircraft as a new concept for European air fleet in 2020. They discussed the design parameters such as maximum stress, ailerons performance, static aeroelasticity and flutter for optimization of the structural weight. Divoux and Frediani 43 studied the flutter characteristics of a lifting system of a PrandtlPlane. They determined the flutter boundary of PrP250 using MSC NASTRAN software. Two different strategies were adopted to increase the flutter boundary; these include increased skin thickness and addition of tip tanks. The results showed that adding 20% thickness to the front wing and 10% to the rear wing, changes the flutter speed to the maximum value. Some new findings on the aeroelasticity of bow wings and PrandtlPlane configurations (a subgroup of joined-wings) are presented by Bombardieri et al. 44,45 Their work shows two contributions: the first case is a critical literature review of the efforts tackling aeroelasticity of box wings and the second contribution is concerned with the aeroelastic analysis considering the antisymmetric modes. Furthermore, they have described the parameters affecting the flutter of box wings such as the location of fuel-tank, rigid-body modes, freeplay of control surfaces, and the front wing tip's region. Silvani and Cavallaro et al. 46,47 investigated the aeroelastic behavior of the PrP250 aircraft. The flutter analysis revealed when composite materials are used in the PrandtlPlane structure, the flutter issues are completely overcome. Recently, Fazelzadeh et al. 48 solved the aeroelastic governing equations of the typical wing section of BWA configuration. They considered plunge and pitch motions for each wing section and modeled the torsion and bending elasticity using two torsional and longitudinal springs; also, the winglet was modeled by a longitudinal spring. The flutter analysis showed that the connection winglet stiffness and angle have significant effect on the box wing flutter speed and frequency. Due to the great importance of the aeroelastic static and dynamic phenomena on the joined-wing configurations, recently, CORIDS (Community Research and Development Information Service) and INEA (Innovation and Networks Executive Agency) have published public information on PARSIFAL project. 49 Therefore, there is not yet a wide literature studying the aeroelasticity of box wings and several aspects need to be understood and analyzed in more depth. In the present study, aeroelastic modeling and flutter analysis of a 3-D box wing aircraft configuration are carried out. The results of this study can help designers to introduce small modifications to the design or come up with new BWA configurations that would avoid the flutter problem.

Governing Equations
A schematic of a BWA configuration is presented in Fig. 1. Because of the geometrical and dynamic complexity of the BWA, several intermediate coordinate systems are applied in this paper. The coordinate system X0Y0Z0 and (XYZ) f, r are fixed on the air plane center of gravity and the root of each wing, respectively. Also the coordinate systems xyz and x'y'z' are located on the root of each un-deformed swept wing and the deformed swept wing, respectively. The transformation between (xyz) f and (xyz) r coordinate systems is as described by Ref. (17). The subscripts f and r denote front and rear wings, respectively.
The box wing tip is shown in Fig. 2. The system consists of two bending and torsion deflections for the front wing (wf, θf) in zf and xf directions as well as two bending and torsion deflections for the rear wing (wr, θr) in zr and xr directions. It should be noted that in this paper, the winglet is modeled by a longitudinal spring and in order to consider the effect of the bending and torsion coupling, two ends of the spring are placed on the gravity centers of the wing sections. It should be mentioned that R and F indexes indicate the center of gravity of the tip section of the rear and front wing before the deformation. A transformation between (xyz)f,r and (x'y'z')f,r coordinate systems is as described by Ref. (19).
in cos sin os sin cos The governing equations of motion are derived via Hamilton's variational principle that can be expressed as Ref. (20).
where Ry and Rz indicate the distance between the center of gravity of the fuselage and the root of wings in y and z directions, respectively. x1, y1 and z1 are the distances of an arbitrary point of the deformed wings from the root. So, the variation of kinetic energy for the left and right wings can be expressed as:   (7) Also, the variation of the strain energy for the left and right wings is as the following equation. (4) 0 where, EI and GJ are bending and torsional rigidity, respectively. To calculate the variation of the spring potential energy, the spring length before (Δr1) and after (Δr2) the deformation should be obtained as The vector of spring deformation and the variation of spring potential energy can be expressed as Eqs. (12) and (13), respectively. 1 2 The variational virtual work of the aerodynamic forces and moments acting on the box wing can be expressed according to the following relations. [ ]  (15) where L and M refer to the aerodynamic lift and moment, respectively. By substituting Eqs. (6)(7)(8)(9) and Eqs. (13)(14)(15) into the Hamilton's variational principle, the PIDEs are derived as Eqs. (16)(17)(18)(19). Note that for every admissible variation (δwf, δθf, δwr, δθr), the coefficient of these variations must be zero.   . .

Wagner Unsteady Model
In order to employ the aerodynamic loadings, Wagner unsteady model is utilized. 24,25 Because of the wing sweep angle, some corrections should be implemented in the aerodynamic model. The following procedure is utilized by Ref. (21, 23).
For the finite-span wings, the modified Strip theory is used to extend 2-D model to 3-D case. Therefore, the following corrections should be applied to the aerodynamic lift and moment. 2 2 .

cos
where A.R. and Clθ are the wing aspect ratio and the lift-curve slope, respectively. The final relations can be expressed after applying the corrections as follows: In these equations, ϕ(t) is Wagner function.

Solution Methodology
Due to the intricacy of the aeroelastic governing equations, the solution is searched by using an approximate solution procedure. Also, to eliminate the time-dependent and parameter-dependent integral terms, some techniques are employed which discussed in this section.

The parameter-dependent integral parts
In order to eliminate the parameter-dependent integral parts, using a new class of generalized functions and by part integral method, the integral terms lead to partialdifferential terms. A new class of generalized functions are derived from hyperbolic tangent family tanh(nx) as follows 50 :

The time-dependent integral parts
In order to remove the time-dependent integral parts, using part integral method and some mathematical simplifications, Eqs. (16)(17)(18)(19) convert into the equations which the timedependent integral terms are eliminated as Refs. 24 and 25.

The numerical approach
Since an analytical solution for the equations of Appendix is difficult, the assumed modes method is employed herein. The bending and torsion deflections (wf, θf, wr, θr) are expanded by means of series of trial functions which only must satisfy geometric boundary conditions. 17 ηi and ψi are time-dependent functions of bending and torsional modes, respectively. Also wi and ϕi indicate mode shapes for wing bending and torsional deflection, respectively, and can be expressed as Ref. 17: ( ) frequency, the maximum and minimum velocities and frequencies are indicated in Table  3.  Table 4. Also, the following non-dimensional parameters are made: 48 where ( ) and are winglet tension stiffness and length, respectively and is the first torsional frequency of the front wing.
A comparison between the approximate calculated mass of the box wing in MSC NASTRAN and Ref. (15) is demonstrated in Table 5. Also, the flutter analyses of the box wing are compared in Table 6. The results show a good agreement between numerical and software results for Ks = 0.5.  The damping variation of the box wing in NASTRAN software, is plotted versus the airstream speed in Fig. 5. Furthermore, in Fig. 6, in order to validate the Wagner aerodynamic model, the swept Goland wing is used. It can be observed that the accuracy of the flutter analysis is desirable and the numerical simulation has good agreement with the published papers.
It is worth noticing that coefficient should be corrected as = 2 / ( �1 + (2 Λ/ ) 2 + 2 Λ) for high sweep angles (Λ > 30°) and considered as = 2 for sweep angles Λ ≤ 30°. In Fig. 7, the effect of winglet tension stiffness on the flutter boundary is exhibited. As it can be observed in Fig. 7a, the stiffness has an optimum value for the flutter speed such that the maximum flutter speed is occurred at K=0.2. Fig.  7b shows that increasing the stiffness increases the flutter frequency, generally. The effect of the box wing sweep angle on the flutter speed and frequency is presented in Fig. 9. The results show that increasing the sweep angle will expand the flutter stability boundary, especially for large sweep angles. But, by increasing the front and rear wings sweep angles, the flutter frequency does not change, significantly. Figure 10 shows the effect of aircraft altitude on the flutter boundary for some selected sweep angles. It can be observed that the flutter speed increases at high altitudes and large sweep angles and the flutter frequency remains almost constant for large sweep angle. Also, effects of the front and rear wing sweep angles on the box wing flutter boundary is demonstrated in Fig. 11 for selected values of the chord ratios. The sweep angle has a minor influence on the flutter speed for the selected values of chord ratio as shown in Fig. 11a. However, increasing the chord ratio enhances the flutter speed and frequency.  Figures 12 and 13 show the effects of bending and torsional rigidity ratio on the flutter boundary for selected values of the winglet stiffness, respectively. According to Fig. 12, increasing the bending rigidity ratio has no effect on the flutter speed at the high winglet stiffness but diminishes the flutter frequency. However, increasing the torsional rigidity ratio significantly decreases the flutter speed and frequency as shown in Fig. 13.
The final parameter that is investigated herein is the box wing span. Fig. 14a demonstrates that adding up the semi-span reduces the flutter speed dramatically at the first and then it remains constant; whereas the flutter frequency shows a completely decreasing trend as shown in Fig. 14b. On the other hand, increasing the winglet tension stiffness decreases the flutter speed and does not change the flutter frequency.

Conclusions
In this paper, mathematical modeling of a 3-D box wing configuration is examined. The resulting aeroelastic partial integro-differential equations are solved by a multi-step procedure and the flutter results are validated with the results of MSC NASTRAN software and previous published papers. The Wagner unsteady model is used to apply the aerodynamic loadings to the equations and in order to analyze the model in MSC NASTRAN, PK method using Doublet-Lattice Subsonic Lifting Surface theory was utilized at the sea-level. Furthermore, the winglet is modeled by a longitudinal spring and the effect of the wing-joint interactions (bending and torsion coupling) is considered by placing two ends of the spring on the gravity centers of the front and rear wing tip sections.
The effects of some design parameters such as winglet tension stiffness, front and rear wings sweep and dihedral angles and aircraft altitude on the flutter boundary are studied.
The analyses revealed that winglet design has a significant influence on the box wing flutter boundary and the following consequences can be resulted: • The maximum flutter speed of the box wing occurs at an optimum winglet stiffness (K=0.2) and on the other hand, a large increase in the winglet stiffness can diminish the flutter velocity. • The dihedral angle of the front wing has no effect on the flutter boundary but the sweep angle of the box wing has a significant influence on the flutter speed, while it has a minor effect on the flutter frequency. • While the chord ratio of the box wing increases, the flutter speed grows up, significantly. • Increasing the box wing span reduces the flutter speed and frequency.