Abstract
Calculating hinge moments during the morphing process is a critical aspect in the folding wing design. The 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. The flat-plate method, which ignores the influence of the airfoil thickness, shows a great simulation error in hinge moment, whereas the high-order panel method can effectively describe the thickness effect and obtain reliable simulation results.
1. Introduction
As a typical morphing scheme, folding wing has received extensive attention in recent years [1–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. The 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.
The 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 mean-plate 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. Thus, 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. Third, this method is insensitive to the paneling scheme; thus, it reduces the requirements for paneling quality, and accurate results are easily obtained. These 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. The influences of airfoil thickness and deficiencies of the flat-plate method are expounded first. Then, 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.
2. Simulation Platform
2.1. Sketch of Folding Wing Aircraft
A typical half-model of a folding wing aircraft is shown in Figure 1. This 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)). The 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.

2.2. Aerodynamic Model
2.2.1. Weakness of the Flat-Plate Aerodynamic Model
The flat-plate aerodynamic model solves the linearized small disturbance equation for subsonic case [11]:where is the freestream Mach number, is the speed of sound, and is the perturbation velocity potential, as defined in the wing frame of reference.
The perturbation velocity potential in equation (1) should satisfy the following wall boundary condition:where is the surface of the aircraft and is the unit out-normal 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 to the x-y plane (take a thin wing that sits near the x-y plane, for instance):where is for the upper surface and is for the lower surface as shown in Figure 2. The upper and lower surfaces can also be expressed by using the mean-plate and the wing thickness , such that

By substituting equation (4) into equation (3), we can obtain the wall boundary condition expressed by and :
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 and a lifting problem , as shown in Figure 3.

The 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:
That 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. Then, 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).

The problem to be solved iswith the wall boundary condition:
The flat-plate method solves this problem by arranging the basic singularity vortex/doublet on the mean-plate and educes the problem of solving the 3D flow domain to a surface solving problem. The 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. The pressure distribution of a folding wing aircraft with an NACA0006 symmetrical airfoil at an angle of attack of 0° and 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. The 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. This 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. Therefore, the majority of the current relevant studies applies theoretical aerodynamic models.
2.2.2. 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 and the unsteady potential :
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 and is of first order for . This is to say that airfoil thickness mainly affects the steady aerodynamic loads, whereas the unsteady loads are unaffected.
Therefore, we can assume that the unsteady part calculated by the flat-plate method is suitable for the flight-folding process and only the steady part 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 established by the high-order panel method and the unsteady part constructed using the flat-plate method is adopted.
2.2.3. 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 , , and . Then, 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. The perturbation velocity potential at the flow field point can be obtained in terms of the source and doublet singularities over the surface and wake , as follows:where and are the source and doublet singularity distribution on surface and wake , respectively (Figure 6). is the distance between the source/doublet and point .

The strengths of and in equation (15) can be determined according to the following boundary conditions:where is the Neumann boundary condition imposed on the exterior surface , is the Dirichlet boundary condition imposed on the interior surface , and is the zero-force condition that should be satisfied by the wake .
In the flat-plate method, the surface is replaced by a mean-plate , and only the doublet singularity is arranged on the mean-plate. The corresponding integral equation is reduced toThe corresponding boundary condition is
Thus, 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). The largest difference from the flat-plate method is that the higher-order panel method is based on the actual object surface. The arranged singularities include doublet that is used to generate lift and source to simulate the thickness effect. Thus, 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. Three sections are selected in this comparison, and the cut plane is y/b = 0.2. The Mach number is 0.3, and the angle of attack is 0°.

(a)

(b)

(c)

(d)
The pressure coefficient of the upper surface at section A is smaller than that of the lower surface, and it is reversed at section C. This finding reflects the influence of the two low-pressure areas on the aerodynamic loading distribution of the fuselage and the outer wings. The pressure coefficient of the upper surface at section B is smaller than that of the lower surface. This 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. Thus, 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. The 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. The aerodynamic modeling process is detailed in Section 2.3.
2.3. Coupling Simulation Platform
The 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. The 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 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. The 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. The flowchart is shown in Figure 10. First, several intermediate configurations are selected. The high-order panel method is used to generate the steady AIC matrices and of the aerodynamic nodes. The doublet lattice method is used to generate the unsteady AIC matrices , , ..., (see Appendix for details). Then, the spline matrix is used to achieve the interpolation between the aerodynamic and structural nodes [14], and the AIC matrices of the structural nodes , , , , ..., are obtained. Finally, the complete AIC matrices are formed by assembling the steady AIC matrices and generated by the high-order panel method and the unsteady AIC matrices , , ..., generated by the doublet lattice method. The AIC matrices of any folding angle are obtained by interpolating the intermediate configurations. The aerodynamic load of the folding wing can be expressed aswhere and are the steady AIC matrices generated by the high-order panel method, , , ..., are the unsteady terms of the AIC matrices generated by the doublet lattice method, is the freestream dynamic pressure, is the reference length, is the freestream velocity, and are the displacement and the aerodynamic force vectors for the structural nodes, respectively, and is the state vector related to and has the following relationship:where 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. The coupling simulation flowchart is shown in Figure 11. The 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. The 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. The 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.
3. Numerical Examples
3.1. Simulation Model
The geometry and dimensions of the folding wing aircraft used in this study are shown in Figure 12. The 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. The 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. The main geometric parameters of the folding wing are listed in Table 1. For more details about the model, refer to [8].

(a)

(b)
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. Thus, 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. The surface network of the folding wing at a folding angle of 120°, which is composed of 2112 panels, is shown in Figure 13(a). The other set is meshed on the mean-plate for the doublet lattice method. The network of the 120° folding angle configuration, which comprises 152 panels, is shown in Figure 13(b).

(a)

(b)
3.2. 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. The 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. The simulation results of the main dynamic parameters are shown in Figure 14.

(a)

(b)

(c)

(d)

(e)

(f)
The results show that the aircraft reaches the initial equilibrium state at approximately 5 s. The initial trimmed angle of attack is 1.44°, and the deflection angle of the aileron is −0.15°. The 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°.
Through 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. The 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.
3.3. Comparison of the Simulation Results of Different Aerodynamic Models
The comparison of the simulation results between the proposed and flat-plate methods [8] is shown in Figure 15. The 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.

(a)

(b)

(c)

(d)

(e)

(f)
The changes in the hinge moments with the folding angle are shown in Figure 16. The 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. The hinge moment reaches a maximum simulation error of 31% when the wings are folded in place. The 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. This simulation error increases with the folding of the wings and reaches a maximum value of 41% when folded in place.

(a)

(b)
The simulation errors caused by the flat-plate method can be explained as follows. The airfoil thickness can cause two low-pressure areas (I and II) between the fuselage and the wings (Figure 17). The 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.

The two low-pressure areas can cause additional normal aerodynamic loads F1 and F2 on the inner and outer wings. F2 (downward) is opposite to the main aerodynamic load (upward) and reduces the hinge moment of the inner and outer wings. F1 (upward) is in the same direction as the main aerodynamic load and increases the hinge moment of the inner wing.
The hinge moment of the inner wing is affected by F1 and F2. When the folding angle is small, the force arm of F2 to the inner wing’s hinge is considerably larger than that of F1. At this time, the influence of F2 dominates despite the larger magnitude of F1 and reduces the hinge moment of the inner wing. When the folding angle increases to a certain level, the force arm of F2 to the inner wing’s hinge becomes as small as that of F1, and F1 is considerably larger than F2. At this time, the influence of F1 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.
The hinge moment of the outer wing is mainly affected by F2. With the folding of the wings, the strength of the low-pressure area II continues to increase, and F2 also increases. Therefore, the hinge moment of the outer wing calculated by the flat-plate method is always larger than that calculated by the proposed method. The 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. The additional aerodynamic loads caused by the two low-pressure 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.
The hinge moments calculated by the two aerodynamic models are compared. The 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.
4. 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)The 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. The results show that the method can effectively consider the thickness effect and obtain reliable calculation results of hinge moments.
Appendix
A. Derivation of the Unsteady AIC Matrices
The 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 is the freestream dynamic pressure, is the AIC matrix, is the Mach number, is the reduced frequency, and and are the aerodynamic force and the displacement vectors, respectively (for the folding wing aircraft containing components in both directions of y- and z-axes).
The AIC matrix 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]. The fitting formula iswhere is the Laplace variable, is the reference length, is the air speed, , , ..., are undetermined matrices, and is the undetermined coefficients. The values in the range of reduced frequency ranges of interest are taken as and is seven for a general fitting accuracy requirement.
In the fitting formula, any mnth element of can be expressed aswhere , , ..., are the seven undetermined coefficients of . Substituting the harmonic condition into equation (A.3), we obtain
The real part of isand the imaginary part is
Assume that , , ..., are the computed AIC matrices using the doublet lattice method with reduced frequencies , , ..., while , , ..., are the mnth elements of , , ..., . The task is to find the free approximation coefficients , , ..., that minimize the total least-square approximation error,where and are the real and imaginary parts of , respectively.
The steady term can be obtained according to the condition ,and the unsteady terms , , ..., are determined using a least-square technique. The result can be expressed aswhere
Figure 18 shows the fitted curves of the aerodynamic influence coefficients of panel A on panel B at and (Figure 13 shows the specific positions of panels A and B). Subscripts y and z represent the components in the y- and z-directions.

(a)

(b)

(c)

(d)
The following error evaluation function is introduced as
Figure 19 shows the fitting errors of matrix, and a maximum percent error is achieved at , . Figure 20 shows the change of 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,
With the introduction of the following state variable:
Equation (A.12) can be expressed as
By transforming equation (A.14) into the time domain, we obtain the final expression of aerodynamic force aswhere is the needed unsteady part of the aerodynamic model. is the steady part and to be replaced by the high-order panel model.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the Fundamental Research Funds for the Central Universities (No. 3082018NS2018060).