Equivalent dynamic modeling of flexible morphing aircraft

Abstract

A morphing aircraft can maintain optimal flight performance by adaptively changing shape. However, large deformation and fast motion of aircraft modules lead to complicated dynamics during morphing process. This paper proposes a mechanistic equivalent model with parameters identified by optimization method. Based on the dynamics of the aircraft modules with large deformation and fast motion and considering the coupling characteristics of rigid bodies and gimbal joints, an equivalent dynamic model of morphing aircraft is built in this study. Considering the huge amount of highly coupled parameters in the equivalent model, particle swarm optimization algorithm is used to identify the equivalent parameters based on the sample data of flexible model. By comparing the simulation results of proposed model to those of rigid model and flexible model, it can be seen that the accuracy of the proposed equivalent model is comparable to that of the flexible model, but the computational load is only 10% of that of flexible model. Further, based on this high-fidelity model with low computational load, an optimized morphing process is obtained, and the attitude variation during morphing is reduced by 4.23%.

Introduction

Morphing aircrafts can achieve optimal flight performance in different tasks and environments by adaptively changing shape. ¹ Through decades of development, scholars have invented many shape-changing modes of morphing aircraft such as variable sweepback, variable wingspan, variable dihedrals, and variable curvature. ² The flight performance of aircraft can be enhanced through morphing, ³ but dynamic characteristics and attitude stability are greatly influenced during morphing process. ⁴ Regarding overall design for a morphing aircraft, the morphing aircraft’s dynamics with different morphing states need to be analyzed. ⁵ In control system design, an accurate model that can sufficiently describe the dynamic characteristics of a morphing aircraft is needed. Modeling of morphing aircraft has become a hot research field. ⁶

The rigid body modeling method has been widely used in morphing aircraft dynamic analysis because of its conciseness. Grant and Lind ⁷ studied time-varying inertias’ effects on an asymmetric variable-sweep morphing aircraft. Tong and Ji ⁸ built a multi-rigid-body dynamic model of an asymmetric variable sweep morphing aircraft via Kane method. Li et al. ⁹ simplified a variable-span and variable-sweep aircraft as a multi-rigid-body system and built the dynamic model by Kane method.

With the use of flexible wings, the structure deformation cannot be ignored, the rigid body models cannot meet the requirement of accuracy, and scholars tend to build flexible multibody model for morphing aircraft. Phani et al. ¹⁰ modeled the morphing aircraft via buckling of structures in multi-element structural components. Snyder et al. ¹¹ investigated the vibration and flutter instabilities during the wing’s folding by finite element method. However, the flexible model’s computation load is too large to analyze flight performance in full flight envelope. For controller design, the model form is too tedious. An accurate modeling method with low computational load is necessary for morphing aircraft.

The equivalent method provides a new insight into dynamic modeling. To relieve the computation load of flexible multibody model, scholars have made attempts to model flexible body with discretized rigid bodies, and such equivalent multi-body modeling method has been rapidly developing. Liang et al. ¹² developed a rigid-flexible coupling model based on the lumped parameter method and analyzed the influence of body number on model accuracy. Chen et al. ¹³ developed a general discretization-based approach to study large deflection of spatial flexible links in compliant mechanisms. The equivalent method shows high accuracy and involves low computational load. For morphing aircraft, equivalent model may relieve computational load and provide a sufficient model form. As a new field, the deformation mechanism shall be explored to build suitable model form for morphing aircraft. In addition, determination of model parameters is also an important issue to be studied. For flexible morphing aircraft, an accurate flexible model has large computation load and complex nonlinear model form, so it is hard to obtain equivalent parameters by experience formulas or analytical formulas.

In this paper, considering the large deformation and fast motion of flexible morphing aircraft, a mechanistic equivalent model composed of rigid bodies and gimbal joints is built. The model can accurately simulate the aircraft’s flexible dynamics with low computational load. Considering the huge amount of coupled parameters of equivalent model, particle swarm optimization method is adopted to obtain the equivalent parameters from the sample data without the need of analytical formulas. After verifying the efficiency and accuracy of the proposed method by simulation, it is applied to morphing process optimization.

The contributions of this paper are listed below:

1. For a sweep-changing aircraft, a novel mechanistic equivalent modeling method composed of rigid bodies and gimbal joints is proposed. The proposed equivalent model has comparable accuracy to that of the flexible model, but the computational load is much less (in fact, the computation load of proposed method is similar to that of the rigid model).

2. For such a model with numerous highly coupled parameters proposed in this work, the model parameters are identified by optimizing the error of equivalent model relative to flexible dynamic model using particle swarm method, which may improve model accuracy compared to empirical or analytical formulas.

3. This equivalent model is applied to morphing phase optimization for the first time. The proposed equivalent model has higher accuracy compared to rigid dynamic model, so the optimized phase can perform better when conducting high-fidelity verification using flexible dynamic model, and the optimization time is acceptable as well. The proposed model is suitable for flight phase design.

Structure of the equivalent model

During the morphing process, the wings rotate around the hinge, and the main deformation forms include bending deformation and torsion deformation. Fast motion between the modules and large elastic deformation of the modules are highly coupled, and the structural dynamics of morphing aircraft is constantly changing. To reflect such deformation phenomenon, each module is equivalent to several rigid bodies that rotate around the spherical sub-axis. In the system dynamic equations, the generalized mass, generalized stiffness and damping matrix vary with the joints’ angle and velocity, which can reflect the time-variant structural dynamics of the morphing aircraft.

Generally, for a flexible body, the bending stiffness and torsion stiffness are different. To embody this feature in the equivalent model, the equivalent stiffness of the three successive rotations in a joint should be different. A gimbal joint can be considered as three successive revolute joints with different stiffness and damping ratio, so it is adopted as the joints between the equivalent bodies within a module. The system dynamic equations are derived by Largrange principle. ¹⁴

In the sweep changing process, each of the wing is equivalent to $l$ rigid bodies connected by gimbal joints. The wing’s rotation along the hinge in the body also provides a degree of freedom. Then in the equivalent model, the degree of freedom of a wing is $n=3(l-1)+1$ , as shown in Figure 1. With the increase of number of equivalent bodies, the accuracy of equivalent model will be improved, while the computational load will increase at the mean time. A suitable number of craft bodies shall be determined considering the efficiency and accuracy of model. We may set $l=3$ initially, and increase $l$ sequentially until the decrement of the model’s error is less than 10% or the increment of the average computational load is more than 10%.

Figure 1.

The equivalent multibody model of a sweep morphing aircraft.

The generalized coordinates are: $\mathit{q}=[{{\mathit{r}}_b}^T,{{\mathit{\Phi }}_b}^T,{\mathit{\lambda }}_1,\cdots ,{\mathit{\lambda }}_{3(l-1)}]$ , where ${\mathit{r}}_b$ is the position of mass center of aircraft body, ${\mathit{\Phi }}_b=[e_1,e_2,e_3]$ is the Euler angle of aircraft body with the rotational sequence of $X,Y,Z$ , ${\lambda }_i$ is the rotation angle of the joint.

The system Lagrangian is:

$$\begin{array}{l}T=\frac{1}{2}{m}_b{\stackrel{.}{r}}_b^T\cdot {\stackrel{.}{r}}_b+\frac{1}{2}{\stackrel{.}{\mathit{\Phi }}}_b^T\mathit{S}\cdot {\mathit{I}}_b\cdot \mathit{S}\cdot {\stackrel{.}{\mathit{\Phi }}}_b+{\sum }_{i=1}^{2l}\left(\frac{1}{2}{m}_i{\stackrel{.}{\mathit{r}}}_i^T\cdot {\stackrel{.}{\mathit{r}}}_i+\frac{1}{2}{\mathit{\omega }}_i\cdot {\mathit{I}}_i\cdot {\mathit{\omega }}_i\right)\\ =\frac{1}{2}{\stackrel{.}{\mathit{q}}}^T\mathit{H}\stackrel{.}{\mathit{q}}\end{array}$$

(1)

where $m_b$ is the mass of aircraft body, ${\mathit{I}}_b$ is the inertia of aircraft body, $m_i$ is the mass of the $i$ th body, ${\mathit{I}}_i$ is inertia of the $i$ th body, ${\mathit{r}}_i,{\mathbf{\omega }}_i$ are the position and angular velocity of the $i$ th equivalent body, which can be obtained from the generalized coordinates and Jacobi matrix.

Based on Lagrange principle, the system dynamics equation is ¹⁵ :

$$\mathit{H}\stackrel{··}{\mathit{q}}+\stackrel{·}{\mathit{H}}\stackrel{·}{\mathit{q}}-\frac{1}{2}{\stackrel{·}{\mathit{q}}}^T\frac{\partial \mathit{H}}{\partial \mathit{q}}\stackrel{·}{\mathit{q}}=\mathit{Q}$$

(2)

where $\mathit{Q}$ is the generalized external force.

To reduce the number of parameters to be identified, each wing is divided into $l$ . bodies in the wing span direction, as shown in Figure 2. $b_1,b_2,\cdots b_{l-1}$ are module division parameters, and $b\in [0,1]$ .The When the parameters of the aircraft material are confirmed, the mass and inertia of each equivalent body are determined by specific division of the module. The equivalent divisions of two wings are assumed to be exactly symmetrical.

Figure 2.

The body division of the equivalent model.

For the morphing aircraft, aerodynamic force and gravity force are main external forces. As for the calculation of aerodynamic force, vortex lattice method has straightforward theoretical explanation and low computational load. ¹⁶ The complex flow field dynamics such as viscosity and hystersis can be represented by additional correction terms with undetermined coefficients. ¹⁷

The generalized force caused by the gravity of aircraft body is:

$${\mathit{Q}}_{\mathrm{gravity}}=m_{b}g{[0,0,1,0\cdots 0]}^T+\sum _{i=1}^n{m}_{i}g{\mathit{e}}_z^T·\frac{\partial {\mathit{r}}_i}{\partial \mathit{q}}$$

(3)

Then the system dynamics equation can be obtained as:

$$\mathit{H}\stackrel{··}{\mathit{q}}+(\stackrel{·}{\mathit{H}}-\frac{1}{2}{\stackrel{·}{\mathit{q}}}^T\frac{\partial \mathit{H}}{\partial \mathit{q}}+\mathit{C})\stackrel{·}{\mathit{q}}+\mathit{Kq}={\mathit{Q}}_{\mathrm{aero}}+{\mathit{Q}}_{\mathrm{gravity}}$$

(4)

where constant diagonal matrix $\mathit{C}$ represents the damping coefficients of the gimbal joints, constant diagonal matrix $\mathit{K}$ represents the stiffness coefficients of the gimbal joints. The proposed model’s equation has simple form and can be integrated to simulate the aircraft system dynamics without the need of finite element analysis, so the proposed model can efficiently reduce the computational load compared with flexible dynamic model.

Parameter identification of the equivalent model

For the proposed model, the undetermined parameters include module division parameter $b_1,b_2,\cdots b_{l-1}$ ( $l-1$ parameters), diagonal matrix of stiffness coefficients $\mathit{K}$ ( $3(l-1)$ parameters) and damping coefficients $\mathit{C}$ ( $3(l-1)$ parameters). There are $7(l-1)$ parameters to be identified. These parameters are highly coupled and hard to be obtained by empirical formula s or analytical formulas. Therefore, the process of parameter identification is carried out. To obtain the sample data, the aircraft dynamics with several sweep changing velocities is calculated by the flexible model and the attitude curve of aircraft body relative to the initial stage is obtained. The attitude change caused by aircraft morphing is important for controller design or morphing phase optimization, so the identification shall mainly focus on the attitude deviation of proposed model from that of flexible model.

The parameters of equivalent model are initialized randomly. Then the proposed equivalent model is evaluated to obtain the attitude curve of the aircraft body. After optimization, the parameters are adjusted to minimize the attitude error of equivalent model according to the identification criterion. For each flight course, the attitude data at time points is obtained, and the attitude deviation can be calculated based on the results of equivalent model in Euler angle form:

$${\mathit{e}}_j\left(i\right)={\left[\begin{array}{ccc}{e}_1-e_{1,\mathrm{sample}}\hfill & e_2-e_{2,\mathrm{sample}}\hfill & e_3-e_{3,\mathrm{sample}}\hfill \end{array}\right]}^T$$

(5)

The identification criterion is to minimize the sum of square of error of $m$ different sweep changing velocities during the whole course:

$$J_{\mathrm{identify}}=\sum _{j=1}^m\sum _{i=1}^n{\Vert {\mathit{e}}_m\left(i\right)\Vert }^2$$

(6)

The optimization method

The equivalent model is highly nonlinear, and the parameters such as dividing length, stiffness, and damping coefficients have different physical properties and orders of magnitude. The traditional gradient optimization method may be trapped in a local minimum or likely to diverge. In this paper, the particle swarm optimization (PSO) method is adopted to operate global search and obtain the optimal model parameters.

As a biology-based evolutionary algorithm, PSO algorithm is derived from the predation behavior of birds. ¹⁸ The algorithm searches the optimal solution from a group of points based on group iterations.

According to the identification criterion, PSO’s fitness function is to minimize the model error $J_{\mathrm{identify}}$ . PSO’s solution space is set as the $7(l-1)$ model parameters. In the optimization process, each particle is randomly initialized with virtual velocity and virtual position in the solution space and modified based on its historical location and current best location. ¹⁹ At each iteration, once individual extreme point and global extreme point, represented as $p_i$ and $p_g$ , are found, the virtual velocity and virtual position of particles can be updated as follows ²⁰ :

$$\begin{array}{l}{v}_{id}^{k+1}=v_{id}^k+C_1{r}_1\left(p_{id}-x_{id}^k\right)+C_2{r}_2\left(p_{gd}-x_{id}^k\right)\\ v_{id}^{k+1}=\{\begin{array}{c}v\\ d_{ma{x}_{id}^{k+}{1}_d^{max}}^{-v_d^{ma{x}_{id}^{k+}{1}_d^{max}}\{}\\ x_{id}^{k+1}=x_{id}^k+v_{id}^{k+1}\end{array}\end{array}$$

(7)

where $v_{\mathrm{id}}^k$ and $x_{\mathrm{id}}^k$ are the current virtual velocity and virtual position of the d th dimension for particle i in the k th iteration, respectively; $C_1$ and $C_2$ represent learning factors; $r_1,r_2$ are random number at the distribution of $(0,1)$ . In the next section, the efficiency and feasibility of proposed method shall be verified, then it will be applied to morphing phase optimization.

Results and discussion

In the simulation, firstly, the attitude curve of the aircraft under three sampled conditions is calculated by the flexible model and the equivalent model’s parameters are identified based on the three samples. Then the proposed model’s efficiency and accuracy are verified under different conditions of constant or uneven sweep changing velocity. Finally, the verified model is applied to design morphing process with less attitude variation. The simulation work flow is shown in Figure 3.

Figure 3.

The workflow of the simulation.

Simulation setting

The morphing aircraft body is simplified as a rigid body with size of 5000 mm (length)  × 700 mm (width)  × 500 mm (height). The size of each wing is 3000 mm (length) × 300 mm(width) × 100 mm (height). The material density of morphing aircraft is set as $\rho =2700\mathrm{kg}/m^3$ . For the flexible model, each wing is divided into five solid beam elements with modulus of elasticity E = 200 MPa, the Poisson’s ratio $\nu =0.3$ ; the stiffness matrix coefficients $\beta ={10}^{-4}s$ . In this paper, the unit of angle is radian. In the morphing process, the flight condition of the morphing aircraft is stayed at velocity = [100,0,0] m/s according to the inertial coordinate system, flight height = 100 m. Initially, the aircraft’s three attitude angles according to the inertial coordinate system are 0 rad. The simulation of the flexible model is carried out in ANSYS software with Intel Core i7 9750H and RAM of 16G.

Parameter identification

In order to obtain the model parameters, the number of samples is set to $m=3$ . The attitude curve of the aircraft is calculated by the flexible model at the sweep changing velocity of 0.2, 0.5, 0.8 rad/s. The time of whole course is set to 2 s. The number of timepoints is set to $n=100$ . Then PSO algorithm is applied to identify the parameters of equivalent model.

The number of craft bodies has a great influence on the accuracy and efficiency of the equivalent model. The parameters of model with different number of craft bodies are optimized, and the results are shown in Table 1. For PSO calculation, the number of particles is set as 2500, the values of both learning factors are set to 1, module division parameters $b$ are restricted by $[0,1]$ , stiffness coefficients $k$ are restricted by $[1,40000]$ and damping coefficients $c$ are restricted by $[1,1000]$ .

Table 1.

The efficiency and accuracy of model with different number of craft bodies.

Number of bodies/ $l$	Identification iteration	Fitness function	Average computational load/s
4	12	$1.56\times 10^{-4}$	0.87060
5	12	$1.329\times 10^{-4}$	0.9495
6	14	$7.37\times 10^{-5}$	1.0125
7	15	$7.048\times 10^{-5}$	1.1617

The table indicates that with the increase of number of equivalent bodies, the accuracy of equivalent model is improved with the sample sweep change velocity. The increase of accuracy between $l=6$ and $l=7$ is less than 5%, while the increase of computational load is more than 10%. Considering the efficiency and accuracy of model, the number of equivalent bodies in this research is set to 6.

The optimized results are used as equivalent model parameters. The attitude curves of the aircraft with the sampled sweep change velocity are obtained and compared to the results of rigid model and flexible model in Figures 4 to 6. The results of rigid model show significant deviation from the results of flexible model. Through contrast, the results of equivalent model with optimized parameters are exactly close to those of the flexible model. This indicates that the flexible characteristics of morphing aircraft shall not be ignored. The proposed equivalent model can show the flexible characteristics in high fidelity. The equivalent PSO method can be used to obtain the parameters suitable for the samples. The performance of the proposed model will be verified in detail later.

Figure 4.

The Euler angle curves of three methods at the sweep change velocity of 0.2 rad/s.

Figure 5.

The Euler angle curves of three methods at the sweep change velocity of 0.5 rad/s.

Figure 6.

The Euler angle curves of three methods at the sweep change velocity of 0.8 rad/s.

Verification of the equivalent model

The equivalent model with optimized parameters shows high accuracy for previous identification of samples. Further, the accuracy of such equivalent model shall be verified in different situations.

Firstly, the sweep changing velocity is set to the interpolation and extrapolation points of the samples: 0.4, 1 rad/s, the attitude curves of the three methods are shown in Figures 7 and 8. The sum of square of attitude error of the whole course of the equivalent model is $1.9523\times 10^{-8}$ with the sweep changing velocity of 0.4 rad/s and $1.0172\times {10}^{-3}$ with the sweep changing velocity of 1 rad/s, indicating high accuracy of equivalent model in different situation.

Figure 7.

The Euler angle curves of three methods at the sweep change velocity of 0.4 rad/s.

Figure 8.

The Euler angle curves of three methods at the sweep change velocity of 1 rad/s.

Then the accuracy of equivalent model is verified under the condition of uneven sweep changing velocity. During the morphing course, the sweeping angle is set as a quadratic function of time:

$$\xi \left(t\right)=a{t}^2+\mathrm{bt}$$

(8)

To obtain the sweeping angle of 1 rad in 2 s, the parameter is set as $a=-0.25,b=1$ . For this situation, the attitude curves of the three methods are shown in Figure 9. With the use of damped elastic joints, the proposed model can also show nonlinearity as the flexible dynamic model. The sum of square of attitude error of the whole course of the equivalent model is $9.0532\times {10}^{-5}$ . The proposed equivalent model shows high accuracy under the condition of uneven sweep changing velocity.

Figure 9.

The Euler angle curves of three methods with uneven sweep change velocity.

For the above two situations, the computation load of three methods is shown in Table 2. While the equivalent model has great accuracy, it’s computation load is just about 10% of the flexible model. The proposed method integrates the high accuracy of flexible model and the low computational load of rigid model.

Table 2.

The computation load of the three methods.

Method	Computational cost with constant sweep change velocity/s	Computational cost with uneven sweep change velocity/s
Flexible model	10.0496	12.8644
Rigid model	0.7432	0.8268
Equivalent model	1.0740	1.1307

Application: Optimize morphing process to reduce attitude variation

The attitude variation caused by the morphing action is seriously concerned. ²¹ After the morphing time and final sweep angle are confirmed, the morphing process shall be optimized to reduce the attitude variation of aircraft body. The equivalent model can simulate the flexible dynamics with low computational load, so it is suitable for calculating the whole-course attitude variation corresponding to each particle’ s parameters in PSO process.

The morphing time is set to 2 s, and the final sweep angle is set to 1 rad. Assume the sweeping angle is a cubic spline function. The interpolation points are set as:

$$\xi \left(0\right)=0,\xi \left(t_1\right)=0.4\mathrm{rad},\xi \left(t_2\right)=0.8\mathrm{rad},\xi \left(2\right)=1\mathrm{rad}$$

(9)

The parameters to be optimized are the time points $t_1,t_2$ . The optimization goal is to minimize the sum of square of attitude variation at the 100 timepoints during the whole course:

$$J_{\mathrm{optim}}=\sum _{i=1}^{100}{\Vert \mathit{e}\left(i\right)\Vert }^2$$

(10)

The parameters are optimized by PSO. The number of particles is set to 50, and the values of both learning factors are set to 1. The fitness function is calculated with rigid dynamic model and the proposed equivalent model, respectively. The optimized morphing phase is shown in Figure 10. The morphing phase optimized with rigid dynamic model, the morphing phase optimized with the proposed equivalent model, and morphing phase under constant sweep changing velocity of 0.5 rad/s are evaluated in the flexible dynamic model. The attitude curves of the three morphing phases are shown in Figure 11.

Figure 10.

The optimized morphing phase.

Figure 11.

Comparation of the attitude curve.

As shown in Table 3, with the equivalent model, the sum of square of attitude variation is reduced by 4.23 % compared to that under constant sweep changing velocity, while the optimized result based on rigid dynamic model is increased by 4.12 %. The proposed method has higher accuracy, so the optimized result can reduce attitude variation more efficiently compared to the rigid dynamic model, and the optimization time is acceptable. The result shows the proposed model is suitable for flight phase design.

Table 3.

Comparison of the optimization.

Method	Converged iterations	Optimization time/s	Result value calculated by flexible model
Optimized by equivalent model	9	513.06	0.0493
Optimized by rigid model	10	416.32	0.0536
Constant sweep changing	–	–	0.0515

Conclusions

With the use of flexible material, the deformation of morphing aircraft will occur inevitably. The large deformation and fast motion of modules cause highly time-variant structural dynamics. The rigid model cannot meet the need of accuracy and the computational load of flexible model is too high. To this end, the paper proposes an equivalent modeling method for morphing aircraft. With rigid bodies and gimbal joints, the nonlinear dynamics caused by large elastic deformation can be drawn in high fidelity. Based on the attitude data calculated by flexible model, the parameters of equivalent model are identified with particle swarm optimization. By comparing the simulation results of proposed model to those of rigid model and flexible model, it can be seen that the accuracy of the proposed equivalent model is comparable to that of the flexible model, but the computational load is only 10% of that of flexible model. Further, based on this high-fidelity model with low computational load, an optimized morphing process is obtained, and the attitude variation during morphing is reduced by 4.23%. Future work will focus on the analysis of the dynamic characteristics of morphing aircraft with more morphing types and under more complicated environment so as to extend the proposed method to more complicated system; moreover, aerodynamic force calculation methods that work well under subsonic/transonic condition shall also be further researched.