Calculation of Hinge Moments for a Folding Wing Aircraft Based on High-Order Panel Method

Calculating hinge moments during the morphing process is a critical aspect in the folding wing design. +e deficiencies of the traditional flat-plate aerodynamic model in the calculation are expounded in this work, and a flight simulation platform based on a high-order panel method is established. On the basis of the platform, a typical flight-folding process of the aircraft is simulated, and the results of different aerodynamic models are compared. Results show that airfoil thickness has a great influence on the aerodynamic loading distribution of wing surfaces and thus affects the hinge moments during the folding process. +e flat-plate method, which ignores the influence of the airfoil thickness, shows a great simulation error in hinge moment, whereas the highorder panel method can effectively describe the thickness effect and obtain reliable simulation results.


Introduction
As a typical morphing scheme, folding wing has received extensive attention in recent years [1][2][3]. Folding wing changes the configuration by folding and expanding the wings during a flight to fit various missions. As a representative issue on folding wing, the calculation of hinge moments during the morphing process has become a research focus.
Lee and Weisshaar [4] calculated the hinge moments at different folding angles through static trim analysis based on the ZONA6 aerodynamic model in ZAERO and analyzed the influences of Mach number and the position of the center of gravity on the results. Reich et al. [5,6] integrated the flexible multibody dynamics method (ADAMS) and the vortex lattice method (an in-house code) to develop an aeroelastic multibody morphing simulation tool for simulating the morphing process of a folding wing aircraft. On this basis, Scarlett et al. [7] conducted a series of wind tunnel simulations on folding wing and investigated the changes in load paths and hinge moments under specific motions. Xu et al. [8] developed the doublet lattice method and deduced the suitable aerodynamic influence coefficient (AIC) matrices for the morphing process. e flight-folding process of a folding wing aircraft was simulated by coupling this aerodynamic model with the multibody structural model in ADAMS, and the general variations of the dynamic parameters (e.g., angle of attack, altitude, deflection angle of ailerons, and hinge moments) during the process were investigated.
e aforementioned studies all applied flat-plate aerodynamic models, such as the vortex lattice method and doublet lattice method, which simulate the lifting surface by arranging the basic singularity vortex/doublet on the meanplate and do not consider the influence of the airfoil thickness. However, this paper discovered that when the folding angle was large, serious aerodynamic interference will occur between the wing surfaces due to airfoil thickness, which will affect the hinge moments during the folding process. us, the flat-plate method shows deficiencies in the calculation of the folding wing's hinge moments.
In comparison with the flat-plate method, high-order panel method has certain advantages [9,10]. First, paneling is based on the aircraft surface, which can effectively simulate the real boundary of the flow field. Second, a quadratic singularity distribution is adopted in the discretization process, which results in high computational precision. ird, this method is insensitive to the paneling scheme; thus, it reduces the requirements for paneling quality, and accurate results are easily obtained. ese characteristics grant the high-order panel method an advantage in the calculation of multibody interference.
In the present work, a calculation method of hinge moments for folding wing based on the high-order panel method is presented. e influences of airfoil thickness and deficiencies of the flat-plate method are expounded first.
en, the high-order panel method is introduced, and a flight simulation platform based on this method is constructed. Finally, the typical flight-folding process of a folding wing aircraft is simulated, and the results of different aerodynamic models are compared.

Sketch of Folding Wing Aircraft.
A typical half-model of a folding wing aircraft is shown in Figure 1. is half-model can be treated as a flexible multibody structure, and each flexible body is a substructure (i.e., fuselage (I), inner wing (II), outer wing (III), and aileron (IV)). e folding angle is defined as the angle between the fuselage and the inner wing. To guarantee the lift performance of the folding wing, the outer wing remains parallel to the fuselage during the folding process.

Aerodynamic
where M ∞ is the freestream Mach number, a ∞ is the speed of sound, and ϕ is the perturbation velocity potential, as defined in the wing frame of reference. e perturbation velocity potential ϕ in equation (1) should satisfy the following wall boundary condition: where S is the surface of the aircraft and n is the unit outnormal vector. In flat-plate method, the first-order approximation of the boundary condition is adopted and the linearized wall boundary condition is transferred from the wing surface S to the x-y plane (take a thin wing that sits near the x-y plane, for instance): where + is for the upper surface S u and − is for the lower surface S l as shown in Figure 2. e upper and lower surfaces can also be expressed by using the mean-plate S c and the wing thickness S t , such that By substituting equation (4) into equation (3), we can obtain the wall boundary condition expressed by S c and S t : Since the wall boundary condition (equation (5)) as well as the small disturbance equation (equation (1)) are linear, it is possible to divide the small-disturbance flow over a thin wing into a thickness problem ϕ t and a lifting problem ϕ c , as shown in Figure 3. e thickness problem is to solve the flow over a symmetric wing with nonzero thickness at zero angle of attack: with the wall boundary condition: According to the wall boundary condition, we can obtain the following result: at is to say that the thickness-induced flow is symmetric about the wing and does not produce lift or lift moment because no pressure differences exist at the corresponding locations of the upper and lower surfaces.
To simplify the modeling, the flat-plate method ignores the airfoil thickness. en, the small-disturbance flow over a thin wing is reduce to a problem of zero-thickness cambered wing at angle of attack-lifting surface ( Figure 4). e problem to be solved is with the wall boundary condition: e flat-plate method solves this problem by arranging the basic singularity vortex/doublet on the mean-plate S c and educes the problem of solving the 3D flow domain to a surface solving problem. e vortex lattice method, doublet-lattice method, and the recently developed unsteady vortex lattice method all belong to the flat-plate method category [12]. Flat-plate method is often used in the steady or unsteady aerodynamic calculation of conventional fixed-wing aircraft, and it has gained numerous achievements.
A folding wing morphing aircraft should complete the folding and unfolding of its wings in flight. Aerodynamic interference will occur between the wing surfaces due to airfoil thickness when the folding angle is large. e pressure distribution of a folding wing aircraft with an NACA0006 symmetrical airfoil at an angle of attack of 0°a nd a folding angle of 120°(Mach number � 0.3) is shown in Figure 5. Two evident low-pressure areas are formed between the fuselage and the inner wing and between the inner and outer wings due to the aerodynamic interference between the wing surfaces. e reason for this formation can be explained as follows. Two airflow channels exist between the fuselage and the wings and narrow at the middle section due to the airfoil thickness. When the airflow passes through the channels, airflow acceleration occurs and low-pressure areas are formed. At this point, the thickness-induced flow is no longer symmetric about the wing, which considerably affects the aerodynamic distribution of the folding wing aircraft and the hinge moments during the morphing process.
When the flat-plate method is considered, the influence of the airfoil thickness is ignored and the calculated aerodynamic loads of the fuselage and the wings should be zero.
is condition indicates the inadequacy of the traditional flat-plate aerodynamic model in the calculation of the folding wing's hinge moments.
Although the CFD method can provide a highly accurate aerodynamic model, the moving grid technique encounters difficulty in handling the wide range of wing motion, and the overlapping grid method cannot automatically fill the continuously increasing gaps between the wings that are covered by flexible skins in the actual model during the folding process [2]. Several technical difficulties in directly applying the CFD method to the simulation of the folding process are still experienced. erefore, the majority of the current relevant studies applies theoretical aerodynamic models.

Modeling Strategy.
For the linearized small disturbance equation (equation (1)), the total perturbation velocity potential ϕ can be expressed as a superposition of the steady potential ϕ 0 and the unsteady potential ϕ 1 : By substituting equation (11) into equation (1), we can obtain the steady and unsteady linearized small disturbance equations: where the thickness effect of the thin lifting surface is of the order for ϕ 0 and is of first order for ϕ 1 . is is to say that airfoil thickness mainly affects the steady aerodynamic loads, whereas the unsteady loads are unaffected. erefore, we can assume that the unsteady part ϕ 1 calculated by the flat-plate method is suitable for the flightfolding process and only the steady part ϕ 0 should be recalculated. In this work, the high-order panel method is used instead of the flat-plate method to calculate the steady term of aerodynamic model, and a modeling strategy of combining the steady part ϕ 0 established by the high-order panel method and the unsteady part ϕ 1 constructed using the flat-plate method is adopted.

Introduction of the High-Order Panel Method.
Similar to the flat-plate method, high-order panel method considers the following steady linearized small disturbance equation in the subsonic case [9]: Let x � x, y � βy, and z � βz. en, we can transform equation (13) into Laplace's equation, as follows: By applying Green's theory and considering the influence of the wake to satisfy the Kutta condition at the trailing edge, equation (14) can be transformed into an integral equation. e perturbation velocity potential at the flow field point (x, y, z) can be obtained in terms of the source and doublet singularities over the surface S and wake W, as follows: where σ and μ are the source and doublet singularity distribution on surface S and wake W, respectively ( Figure 6). r is the distance between the source/doublet and point (x, y, z). e strengths of σ and μ in equation (15) can be determined according to the following boundary conditions: where (zϕ 0 /zn) � −V ∞ · n is the Neumann boundary condition imposed on the exterior surface S + , ϕ 0 � 0 is the Dirichlet boundary condition imposed on the interior surface S − , and (zϕ 0 /zx) � 0 is the zero-force condition that should be satisfied by the wake W.
In the flat-plate method, the surface S is replaced by a mean-plate S c , and only the doublet singularity is arranged on the mean-plate. e corresponding integral equation is reduced to e corresponding boundary condition is us, the high-order panel method is to solve integral equation (15) under the boundary condition in equation (16), whereas the flat-plate method is to solve the integral equation (17) under the boundary condition in equation (18). e largest difference from the flat-plate method is that the higher-order panel method is based on the actual object surface. e arranged singularities include doublet that is used to generate lift and source to simulate the thickness effect.
us, the high-order panel method theoretically considers the thickness of the airfoil. Furthermore, the higher-order panel method (e.g., the open source software PANAIR used in this work) adopts a quadratic singularity distribution in the discretization process of equation (15), which indicates a high accuracy.
A comparison of the steady pressure distribution calculated by the higher-order panel (PANAIR) and CFD methods (i.e., Euler's result) is shown in Figure 7. ree sections are selected in this comparison, and the cut plane is y/b � 0.2. e Mach number is 0.3, and the angle of attack is 0°.
e pressure coefficient of the upper surface at section A is smaller than that of the lower surface, and it is reversed at section C. is finding reflects the influence of the two lowpressure areas on the aerodynamic loading distribution of the fuselage and the outer wings. e pressure coefficient of the upper surface at section B is smaller than that of the lower surface. is finding indicates a greater influence of the low-pressure area between the fuselage and the inner wing than that between the inner and outer wings. In addition, the calculation results based on the high-order panel method and CFD method agree well. us, the high-order panel method can effectively consider the effect of airfoil thickness.
It should be noted that all existing software of high-order panel method can solve a steady flow field. e boundary condition considered is a steady condition, and the calculated strength of the singularity is a steady strength, which is only applicable to the calculation of a steady flow field and cannot be used directly in the folding process.
According to the modeling strategy in Section 2.2.2, we can do a combination of the high-order panel method and flat-plate method to obtain a complete aerodynamic model that contains the steady and unsteady terms. e aerodynamic modeling process is detailed in Section 2.3.

Coupling Simulation Platform.
e simulation platform in this study is based on study in [8]. Certain necessary modifications are applied to the aerodynamic model to consider the thickness effect. e platform is used to realize the dynamics simulation of the flight-folding process and contains three main modules, namely, flexible multibody dynamic modeling, unsteady aerodynamic load modeling, and flight control modeling ( Figure 8).
In structural modeling, the platform integrates the capabilities of ADAMS and NASTRAN, thereby combining advanced multibody dynamics models and finite element models to simulate complex flexible multibody systems [13]. First, through NASTRAN software, the classic

Mathematical Problems in Engineering
Craig-Bampton modal synthesis method is used to perform a modal analysis for each substructure to generate a modal neutral file suitable for importing into the ADAMS modeling environment. Each substructure is then sequentially introduced into ADAMS and assembled using revolute joints. e assembled full model of a folding wing is shown in Figure 9.
In aerodynamic modeling, the steady aerodynamic model established by the high-order panel method and the unsteady part constructed using the flat-plate method are assembled to obtain a complete aerodynamic model. e flowchart is shown in Figure 10 en, the spline matrix is used to achieve the interpolation between the aerodynamic and structural nodes [14], and the AIC matrices of the structural nodes c i g , E i g , A i,1 g , A i,2 g , ..., A i,7 g are obtained. Finally, the complete AIC matrices are formed by assembling the steady AIC matrices c i g and E i g generated by the high-order panel method and the unsteady AIC matrices A i,2 g , A i,3 g , ..., A i,7 g generated by the doublet lattice method. e AIC matrices of any folding angle θ are obtained by interpolating the intermediate configurations. e aerodynamic load of the folding wing can be expressed as where c θ g and E θ g are the steady AIC matrices generated by the high-order panel method, A θ,2 g , A θ,3 g , ..., A θ,7 g are the unsteady terms of the AIC matrices generated by the doublet lattice method, q ∞ is the freestream dynamic pressure, b is the reference length, V is the freestream velocity, u g and f g are the displacement and the aerodynamic force vectors for the structural nodes, respectively, and y i is the state vector related to u g and has the following relationship:     where k i is a selected reduced frequency in the range of interest.
Based on the capability of secondary development, aerodynamic loading and updating are realized in ADAMS. e coupling simulation flowchart is shown in Figure 11. e subroutine calls the SYSARY and SYSFNC macros to read the position, velocity, and acceleration information of the structural nodes from the ADAMS/ Solver and then passes this information to the aerodynamic calculation program. e aerodynamic load is calculated through the AIC method and then fed back to the ADAMS/Solver.
In addition, the folding process is accompanied by large changes in mass distribution and aerodynamic load distribution. To maintain flight stability during the folding process, we establish a longitudinal stabilization control system in the ADAMS/Controls module. e control system takes the pitch angle and altitude as feedback signals and uses a simple PID control method to drive the ailerons and maintain the longitudinal stability of the aircraft.

Simulation Model.
e geometry and dimensions of the folding wing aircraft used in this study are shown in Figure 12. e fuselage and wings are composed of skin and inner skeleton and assume an NACA0006 airfoil shape. An aileron is set at 80%-100% of the outer-wing chord and 11%-83% of the outer-wing span and serves as a control surface during the flight-folding process. e intersections of the spars and ribs are selected as the interpolation points of the structure, which are used for the interpolation between the aerodynamic and structural nodes. e main geometric parameters of the folding wing are listed in Table 1. For more details about the model, refer to [8].
According to Section 2.3, the complete aerodynamic model of the folded wing is assembled from the steady terms calculated by the high-order panel method and the unsteady terms determined by the doublet lattice method. us, we need to establish two sets of aerodynamic models for the aircraft. One set is meshed on the actual wing surface for the high-order panel method. e surface network of the folding wing at a folding angle of 120°, which is composed of 2112 panels, is shown in Figure 13(a). e other set is meshed on the mean-plate for the doublet lattice method. e network of the 120°folding angle configuration, which comprises 152 panels, is shown in Figure 13(b).

Simulation of Flight-Folding Process.
On the basis of the flight simulation platform constructed in this study, the typical flight-folding process of the folding wing aircraft is simulated, and the general variation of the dynamic parameters during the process is investigated.
In the simulation, the initial flight altitude of the aircraft is set to 2 km, and the Mach number is 0.3. e aircraft maintains an unfolded configuration in the initial 10 s for the preliminary trimming calculation and begins to fold at 10 s, following the cosine law. After 30 s, the aircraft is folded in place; such folded configuration is held to maintain a level flight. e simulation results of the main dynamic parameters are shown in Figure 14.
e results show that the aircraft reaches the initial equilibrium state at approximately 5 s. e initial trimmed angle of attack is 1.44°, and the deflection angle of the aileron is −0.15°. e wings begin to fold at 10 s. As the folding angle increases, the effective span and area of the wing continues to decrease, and the aircraft drops. To maintain a longitudinal balance, the ailerons are controlled to deflect upward, thereby generating a head-up torque that increases the angle of attack to maintain the balance of the lift. At 40 s, the wing is folded in place and quickly reaches a new equilibrium state, which is the trimmed condition of the folded wing configuration. In comparison with the initial equilibrium state, the angle of attack is increased to 5.27°, and the ailerons deflect upward to −4.78°.
rough the simulation, we can monitor the changes in hinge moments during the folding process. In the initial unfolded state, the driver of the inner wing should drive the inner and outer wings to move together, and the hinge moment is relatively large at 2697 N m. e hinge moment of the outer wing in this state is relatively small at 931 N m. During the folding of the wings, the hinge moments of the inner and outer wings tend to initially decrease and then increase. When folded into place, the hinge moments of the inner and outer wings increase to 3391 N m and 1552 N m, respectively.

Comparison of the Simulation Results of Different Aerodynamic Models.
e comparison of the simulation results between the proposed and flat-plate methods [8] is shown in Figure 15. e flat-plate method does not consider the influence of the airfoil thickness, which results in a larger simulation error compared with the proposed method, especially in the calculation results of the hinge moments.  Mathematical Problems in Engineering e changes in the hinge moments with the folding angle are shown in Figure 16. e hinge moment of the inner wing calculated by the flat-plate method is larger than that calculated by the proposed method before the folding angle reaches approximately 90°folding angle and smaller after this folding angle is exceeded. e hinge moment reaches a maximum simulation error of 31% when the wings are folded in place. e hinge moment of the outer wing calculated by the flat-plate method is always larger than that calculated by the proposed method during the entire folding process. is simulation error increases with the folding of the wings and reaches a maximum value of 41% when folded in place. e simulation errors caused by the flat-plate method can be explained as follows. e airfoil thickness can cause two low-pressure areas (I and II) between the fuselage and   the wings (Figure 17). e low-pressure areas have three characteristics. First, the strength of the two low-pressure areas increases with the folding of the wings. Second, the strength of area I is always higher than that of area II. Lastly, when the folding angle is larger than 90°, the left and right inner wings are close to each other, which causes a severe aerodynamic interference and further enhances the strength of low-pressure area I. e two low-pressure areas can cause additional normal aerodynamic loads F 1 and F 2 on the inner and outer wings. F 2 (downward) is opposite to the main aerodynamic load (upward) and reduces the hinge moment of the inner and outer wings. F 1 (upward) is in the same direction as the main aerodynamic load and increases the hinge moment of the inner wing. e hinge moment of the inner wing is affected by F 1 and F 2 . When the folding angle is small, the force arm of F 2 to the inner wing's hinge is considerably larger than that of F 1 . At this time, the influence of F 2 dominates despite the larger magnitude of F 1 and reduces the hinge moment of the inner wing. When the folding angle increases to a certain level, the force arm of F 2 to the inner wing's hinge becomes as small as that of F 1 , and F 1 is considerably larger than F 2 . At this time, the influence of F 1 dominates and increases the hinge moment of the inner wing. Near the end of the folding process, the strength of low-pressure area I is greatly enhanced, which significantly increases the hinge moment of the inner wing.
e hinge moment of the outer wing is mainly affected by F 2 . With the folding of the wings, the strength of the lowpressure area II continues to increase, and F 2 also increases. erefore, the hinge moment of the outer wing calculated by the flat-plate method is always larger than that calculated by the proposed method. e simulation error obtained by the flat-plate method continues to increase with folding and reaches a maximum value when folded in place. Low-pressure area I is located on the upper surface, whereas low-pressure area II is located on the lower surface. e additional aerodynamic loads caused by the two lowpressure areas can offset each other, which reduces the effect on the overall load of the aircraft; such effect is reflected in the small simulation errors in terms of angle of attack, altitude, and deflection angle of the aileron calculated through the flat-plate method. e hinge moments calculated by the two aerodynamic models are compared. e flat-plate method ignores the influence of the airfoil thickness and causes a large simulation error of the hinge moment. In conclusion, the proposed approach based on the high-order panel method can effectively describe the thickness effect and obtain reliable calculation results of hinge moments.

Conclusions
(1) Airfoil thickness has an important influence on the aerodynamic load distribution of the folding wing and on the hinge moments of the folding process. (2) e flat-plate aerodynamic model does not consider the influence of the airfoil thickness, which will cause a large simulation error of the hinge moment. (3) A calculation method of the hinge moments for folding wing based on the high-order panel method is established. e results show that the method can effectively consider the thickness effect and obtain reliable calculation results of hinge moments.

A. Derivation of the Unsteady AIC Matrices
e aerodynamic model for the folding wing with a specific and fixed folding angle can be formulated using the doublet lattice method as follows [15]: where q ∞ is the freestream dynamic pressure, Q is the AIC matrix, M ∞ is the Mach number, k is the reduced frequency, and f and u are the aerodynamic force and the displacement vectors, respectively (for the folding wing aircraft containing components in both directions of y-and z-axes). e AIC matrix Q is obtained in simple harmonic conditions and cannot be directly applied to the aerodynamic calculations in the time domain. To calculate the unsteady aerodynamic force in arbitrary motion, the Roger method is used in this study to derive a rational approximation of the unsteady AIC matrices in the Laplace domain [16,17]. e fitting formula is where s is the Laplace variable, b is the reference length, V is the air speed, A 1 , A 2 , ..., A n are undetermined matrices, and r j−3 is the undetermined coefficients. e n − 3 values in the range of reduced frequency ranges of interest are taken as r j−3 and n is seven for a general fitting accuracy requirement.
(A.11) Figure 19 shows the fitting errors of Q matrix, and a maximum percent error ε max � 0.79% is achieved at m � 2, n � 280. Figure 20 shows the change of ε max with the folding angle. For the entire folding range, the percent errors are within 1%, indicating high fitting accuracy.
Substitution of equation (A.2) into equation (A.1) yields the fitted aerodynamic model, (A.12) With the introduction of the following state variable: A j y j−3 . (A.14) By transforming equation (A.14) into the time domain, we obtain the final expression of aerodynamic force as where A 2 _ u(b/V) + A 3 € u(b 2 /V 2 ) + 7 j�4 A j y j−3 is the needed unsteady part of the aerodynamic model. A 1 u is the steady part and to be replaced by the high-order panel model.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest regarding the publication of this paper.