Aerodynamic, Kinematic and Control Coupling Simulation of Aircraft Maneuvering Flight ()
1. Research Background
The traditional approach to control system design follows an iterative process of “test flight and improvement” [1]. Initially, the unsteady aerodynamics problem and rigid body dynamics problem are decoupled based on the linear small disturbance hypothesis. Through wind tunnel testing, static and dynamic data of the aircraft are acquired to establish an aerodynamic database. Subsequently, an aerodynamic model is developed through parameter identification. Based on this model, the control system design is accomplished with further enhancements and optimizations made during flight testing. However, this design methodology is typically suitable for low angles of attack and quasi-steady flight processes. When applied to high angles of attack and highly maneuverable flights, a significant nonlinear cross-coupling effect emerges between the aircraft’s motion and the flow around it, rendering the linear small disturbance hypothesis inapplicable. Moreover, due to neglecting unsteady aerodynamic characteristics and multidisciplinary coupling issues, ensuring flight safety often requires amplifying uncertainties in aerodynamic data during flight control design or artificially constraining the controllable safety range. Consequently, this results in compromised flight performance and an expensive and time-consuming overall design iteration process.
Currently, the aerodynamic-motion-control coupling solution method has been preliminarily researched and applied in China, primarily concentrated at institutions such as the China Aerodynamics Research and Development Center, China Academy of Aerospace Aerodynamics Technology, and relevant aviation colleges and universities. Among them, CHEN Qi [2] was the pioneer in investigating the pitch motion rule and aerodynamic characteristics of a square-section missile under open-loop control. Additionally, he conducted initial research on integrating numerical simulation technology of unsteady flow field, aircraft motion, and rudder deflection control based on dynamic grid technology. The follow-up research further delved into the aerodynamic control integrated maneuvering flight [3] [4]. In order to simulate missile flight numerically, Da Xingya et al. [5] employed a nested grid and Adams predictive correction method to achieve the coupled calculation of flight mechanics equations and flow field control equations. They developed an aerodynamic/flight mechanics numerical calculation method and software, and subsequently validated and optimized the PID control law for longitudinal balancing flight of missiles based on this foundation. The influence of closed-loop control parameters is also investigated [6]. Chang Xinghua et al. [7] integrated the solution of rigid body dynamics equations and fluid mechanics control equations with flight control laws, establishing a virtual numerical simulation method for aircraft and validating it. Wang Xiaobing [8] conducted numerical simulations to investigate the aerodynamic/kinematic coupling problem in missile longitudinal maneuvers, taking into account the impact of roll instability on missile pitching maneuvers. LAIPING Z, et al. [9] developed a numerical virtual flight solver based on computational fluid dynamics (CFD) technology, validated the accuracy of the numerical virtual flight system through wind tunnel experiments, and successfully applied it to missile control rate design [10]. Additionally, this solver was employed to numerically simulate quasi-cobra maneuvers of fighter aircraft [11]. XI Ke et al. [12] developed a virtual flight numerical simulation platform capable of solving the unsteady N-S equation, rigid body six-degree-of-freedom motion equations, and coupled rudder deviation control law. They conducted virtual flight numerical simulations on the fundamental winged missile shape. Building upon this technology, HUANG Yu et al. [13] simulated the multi-degree-of-freedom unsteady motion, controlled characteristics, and control parameter settings of both the return capsule of “Origin” and the basic winged missile shape.
Currently, the primary focus of virtual flight research both domestically and internationally is on missiles, followed by fighter jets and hypersonic vehicles. However, most of these simulations involve simple shapes, with relatively few complex aircraft shapes being simulated. The main emphasis lies in studying the pitching and rolling motion of aircraft with single degrees of freedom, while other maneuvering actions are less frequently simulated. As a multidisciplinary coupling technology, further development is required for numerical virtual flight technology to effectively simulate highly maneuverable aircraft. Advanced aerodynamic and motion-control integrated coupling methods should be developed to advance aircraft design and performance evaluation, thereby driving the progress of new aviation technologies. The study of aircraft motion and control mechanisms holds significant theoretical importance as well as promising practical applications.
2. Technical Method
In order to address the coupling problem between aerodynamics and kinematics of aircraft during flight, this study aims to establish a virtual flight simulation platform that incorporates the numerical solution method for unsteady airflow, unstructured nested grid technology, and numerical solution method for flight mechanics equations. Additionally, it integrates the control of aircraft rudder surface deflection to achieve an integrated design of motion and control for capturing the aerodynamic characteristics of the aircraft. The key technologies involved in developing the solver and virtual flight platform are as follows:
2.1. CFD Numerical Solution Method
Dynamic grid generation is an important work in the calculation of unsteady flow field with geometric deformation or relative motion of many bodies. With the boundary movement, the flow field calculation grid, including the surface grid and the volume grid, must be updated at each time step to adapt to the change of the flow field calculation domain. Obviously, it is not feasible to re-generate the grid at each step with a large amount of computation. Therefore, dynamic unstructured nested grid technology is adopted in this paper to divide the CFD calculation grid into dynamic unstructured nested grid, and the aircraft flow field is divided into multiple regions with overlapping components. Each region generates an independent grid and is solved separately [14]. The changes on the geometric boundary of the model are reflected on the internal grid by algebraic interpolation, so that the information interaction between different flow fields is established.
In the aspect of CFD numerical solution, three-dimensional unsteady Reynolds means N-S equation (RANS equation), which is currently developed and highly reliable, is used to calculate and simulate the static and dynamic aerodynamic data of aircraft. The unsteady, compressible flow is discretized in two dimensions, space and time. For the space, the finite volume method is used to compute the space discretization. For time, the pseudo-time step is discretized by a one-order backward difference. The physical time step is an implicit second-order backward difference scheme. In addition, in order to make Reynolds mean N-S equation closed, S-A (Spalart-Allmaras) turbulence equation needs to be introduced. The S-A equation model is simple to calculate, and the boundary conditions of inflow and outflow and wall surface are also convenient to implement.
2.2. Flight Dynamics Model
The motion law of aircraft under the action of external forces is generally described by the equation of motion, that is, the form of differential equation is used to describe the change law of the motion and state parameters of aircraft with time. The basic equation of aircraft motion, which can be divided into dynamic equation and kinematic equation, of which six dynamic equations describe the motion of the center of mass and the rotation of the aircraft around the center of mass, and six kinematic equations describe the change of the position and attitude of the aircraft in space. Before establishing the equation of aircraft motion, we must first assume that the aircraft is rigid, and the mass is constant. Its dynamic equation is as follows:
(1)
where,
represents the component of flight speed on the x,y, and z axes respectively;
represents the component of the rotational angular velocity of the aircraft in each coordinate axis of the aircraft body coordinate system, representing the roll Angle velocity, pitch Angle velocity and yaw angular velocity respectively.
represents the projected components acting on the total aerodynamic force of the aircraft and the T-direction moving coordinate system of the engine thrust, where
and
are pitch Angle and roll Angle, respectively.
are the components of the aerodynamic moment of the aircraft on each coordinate axis, respectively representing the rolling moment, pitching moment and yawing moment. Coefficients
in the formula are respectively expressed as:
The kinematics equation of the aircraft is as follows:
(2)
In the formula,
are respectively the rolling angle, pitch angle and yaw angle of the aircraft, and
represents the flight speed. Equations (1) and (2) together determine the nonlinear functional relationship between the state variable
and the control input vector
. All the 12 equations described are closed. Since the Runge-Kutta method is very suitable for solving the initial value problem of differential equations with constant coefficients, the value of the required state vector is the value at the beginning of the time step. Therefore, this paper adopts the fourth-order Runge-Kutta method to solve the parameters such as angular velocity and position at a given time step, and the calculation formula is as follows:
(3)
2.3. Aircraft Control System Design Method
Figure 1 shows the schematic diagram of the closed-loop attitude control system of the aircraft in this paper. In MATLAB simulation, the dynamics model of the aircraft is in the form of a simplified transfer function, while in the numerical virtual flight simulation platform, it is replaced by unsteady flow numerical calculation and numerical solution of the six-degree-of-freedom motion equation. The control method adopts PID control, which is relatively mature and widely used at present, to control the attitude Angle of the aircraft. In this paper, discrete positional PID and incremental PID control algorithms are used to discretize Equation (4). The sampling moment point
is used to represent continuous time, the numerical integral of matrix method is used to replace the integral, and the first-order backward difference is used to replace the differential. The discretized positional PID expression can be obtained:
(4)
![]()
Figure 1. Schematic diagram of a closed-loop attitude control system.
The integral term of positional PID represents the accumulation of errors, and the resulting output value is highly influenced by past states. When the actuator requires incremental control rather than absolute position values, an incremental PID control algorithm can be employed. In comparison to positional PID, incremental PID eliminates error accumulation and reduces the impact of calculation errors on determining the control quantity. The discrete expression for incremental PID is as follows:
(5)
where,
(6)
Since the deflection rate and Angle of the aircraft’s rudder surface are limited, the deflection speed and angle are respectively limited here, and the integral limiting is increased to prevent the integral term of the position-type PID from reaching saturation. The limiting Settings are as follows:
1) Position limiting: the rudder surface of the aircraft can only be deflected within the range of umin~umax, if u > umax, then u = umax; If u < umin, u = umin.
2) Speed limiting: Keep the deflection speed of aircraft rudder surface between rmin and rmax, and the deflection speed is expressed as
. The speed-limiting expression is as follows:
(7)
2.4. CFD, RBD and FCS Integrated Coupling Strategy
In this paper, the loose coupling method is used to realize the coupled solution. The loose coupling method can keep the independent solution of each module and reduce the difficulty of programming. The coupling diagram is shown in Figure 2, that is, the N-S equation, the six-degree-of-freedom motion equation
and the flight control system are solved independently in each time step, and the coupled solution between the modules of various disciplines is realized through the staggered time steps in the time domain.
3. Application of Numerical Virtual Flight Platform in Fighter Aircraft
In this part, the aerodynamic, motion and control numerical calculation methods will be used to carry out simulation. It includes the calculation of static and dynamic aerodynamic characteristics and the design of a nonlinear controller. Using the established virtual flight simulation platform, by adjusting PID controller parameters, the flight state of F16 fighter jet at high angle of attack is simulated. The flight trajectory simulation of multi-channel roll pull-up coupling flight mode of F16 is also carried out.
3.1. Models and Grids
The calculation model adopted in this paper is the F-16 fighter model, as shown in Figure 3.
The calculation domain and grid structure of F-16 fighter model are shown in Figure 4. The simplified fighter model retains the intake port, leading edge slats,
![]()
Figure 4. Calculation domain and grid structure.
ailerons, elevators and rudder, so as to restore the aerodynamic shape of the aircraft model as much as possible. In the process of grid division, the wings, the front and rear edges of the vertical tail and the tail, and the inlet are encryption. Using the elevator as an example, the grid nesting relationship is shown in Figure 5.
![]()
Figure 5. The nesting relationship elevator.
The grid contains six sets of sub-grids, namely, one set of background grid, one set of body grid, two sets of control elevator grid and two sets of control aileron grid, and the total number of grid units is about 20 million.
3.2. Elevator Control of Pitch Channel
In order to realize the pull-up control of the fighter model at large angle of attack, the gain scheduling design method based on feature points is adopted in this paper, and multiple feature points are selected within the controlled angle of attack range. A linear PID controller is designed for each feature point, and these controllers are combined to realize nonlinear control of high angle of attack maneuvering process. In this paper, aiming at the pull-up process of the fighter from 0˚ to 40˚ angle of attack, a total of 8 feature points are selected every 5˚ angle of attack to design PID controllers respectively. The schematic diagram of feature point selection and control system design is shown in Figure 6.
![]()
![]()
Figure 6. Feature point selection and control system design diagram.
On this basis, a numerical virtual flight simulation is carried out for the pull-up process at different initial angles of attack, and the step-up process with an incremental angle of attack of 5˚ is realized, respectively. The maximum time step is 5 ms according to the minimum grid scale, incoming flow velocity and Currant number, so the time step of the following numerical simulation is 1 ms, the maximum elevator deflection angle is 25˚, and the maximum deflection speed is 30˚/s. When the initial angle of attack is 0˚, the variation curve of the aircraft’s angle of attack and rudder yaw is shown in Figure 7. It can be seen that under the action of PID controller, the angle of attack can reach the control target, and the step-up of the incremental angle of attack of 5˚ can be achieved. The adjustment time of the whole process is 0.49 s.
Figure 8 shows the change curves of the longitudinal aerodynamic force, torque and pitch angle velocity of the fighter in the process of pull-up. During the lift-up of the aircraft angle of attack from 0˚ to 5˚, the time is from the initial moment to 0.11s, the elevator of the aircraft is deflecting upward, the elevation moment coefficient increases, the elevation angle velocity increases linearly, the angle of attack of the aircraft is rapidly pulled up, and the lift coefficient and drag coefficient both increase gradually.
![]()
(a) (b)
Figure 7. When the initial angle of attack is 0˚, the change curve of the angle of attack and the rudder angle. (a) Angle of attack; (b) Elevator deflection angle.
Figure 9 shows the distribution of pressure coefficients on the longitudinal symmetric surface and upper and lower surfaces of the aircraft during the process of the angle of attack from 0˚ to 5˚. It can be seen from the figure that during the process of pulling up, the head, cockpit and intake port of the fighter jet exhibit significant high-pressure regions, and due to the outflow of air from the tail nozzle of the fighter jet, a large pressure is also generated at the tail. From the initial moment to 0.11 s, the relative angle between the elevator and incoming flow is negative due to the linear upward deflection speed of the elevator, so the pressure coefficient on the upper surface of the elevator slightly increases, while the pressure coefficient on the lower surface gradually decreases, resulting in the moment of raising the head of the aircraft, the angle of attack of the aircraft is gradually raised, and the negative pressure zone on the upper surface of the wing is significantly increased. During the period from 0.11 s to 0.21 s, the elevator began to tilt down, the aircraft angle of attack continued to increase, the negative pressure area on the upper surface of the elevator expanded, the lower surface pressure coefficient increased, the upper surface pressure coefficient of the edge wing and the wing decreased, and the lower surface pressure coefficient increased. At 0.21 s, the elevator reaches the initial position and then continues to deflect downward. The pressure coefficient on the upper surface of the elevator gradually decreases, the pressure coefficient on the lower surface increases, the angle of attack increases, and the negative pressure area on the upper surface of the aircraft wing and wing gradually expands.
3.3. Elevator Control of Pitch Channel at High Angle of Attack
In the previous section, the same research method used for 0˚ - 5˚ pull-up was adopted, and 8 groups of control results of 5˚ - 10˚, 10˚ - 15˚, 15˚ - 20˚, 20˚ - 25˚, 25˚ - 30˚, 30˚ - 35˚ and 35˚ - 40˚ were obtained successively. The specific control parameters were shown in Table 1. The process of high angle of attack from 0˚ to 40˚ is realized.
Figure 10 shows the variation curves of the angle of attack and the elevator declination of the fighter in the process of pull-up. Figure 11 shows the variation
![]()
Figure 9. Angle of attack 0˚ to 5˚ pull up process flow field changes.
![]()
Table 1. PID controller parameters at each angle of attack.
(a)
(b)
Figure 10. 0˚ - 40˚ angle of attack during pull-up, the angle of attack and rudder angle change. (a) Angle of attack; (b) Elevator deflection angle.
curves of aerodynamic force, moment and pitch angle speed in the process of pull-up.
It can be seen that from the initial moment to 0.426s, the elevator keeps deflecting upward until −15.09˚, the angle of attack is rapidly pulled up, the pitch angle speed increases from 0 to 2.3 rad/s, the lift coefficient and the drag coefficient are gradually increased, the tilt moment coefficient increases slowly, and the roll moment coefficient and side force coefficient change little. However, from 0.426 s to 0.455 s, the rudder of the aircraft has a short downward deflection, which makes the lift and pitching moment provided by the rudder surface decrease, the pulling speed of the angle of attack slow down, the lift coefficient and the drag coefficient slowly increase, and the pitching moment coefficient of the lower head gradually increase. From 0.455 s to 0.6 s, although the elevator continuously deflections upward from −14.52˚ to −18.87˚, the angle of attack of the aircraft changes little. During this process, the pitch moment provided by the elevator can only make the angle of attack of the aircraft slightly raised, and the pitch Angle speed quickly drops to 0˚. The lift coefficient and drag coefficient gradually decrease, and the bow moment coefficient gradually decreases. During this process, the oscillation amplitude of rolling moment coefficient and side force coefficient gradually increases. At 0.59s, the head down moment began to change into the head-up moment, and then from 0.6 s to 0.734 s, the elevator surface continued to deflect upward, the head-up moment coefficient slowly increased, the pitch angle speed gradually increased, and the angle of attack was pulled up again. At 0.734 s, the rudder surface began to deflect downward, and then slowly deflected and finally maintained at −18.24˚. The angle of attack increased slowly and entered the error zone at 0.933 s, and finally remained at 40˚. The lift coefficient and drag coefficient finally oscillated periodically near a fixed value, and the pitch moment coefficient and pitch angle velocity oscillated steadily near 0. The oscillation amplitude of rolling moment coefficient and side force coefficient decreases gradually and oscillates periodically with the stabilization of the angle of attack.
Figure 12 shows the eddy current structure, cross-section vorticity and space
Figure 12. The variation of vortex, cross-section vorticity and space flow line during 0˚ to 40˚ angle of attack.
streamline changes of the aircraft during high angle of attack. It can be seen that obvious vortices appear on the upper surface of the aircraft, including edge wing vortices, wing tip vortices, leading-edge vortices of the main wing, and horizontal tail vortices. With the increase of the angle of attack, the range of vortex structure generated by the edge strip wing of the aircraft becomes larger. At t = 0.3 s, the detached vortex generated by the edge strip wing extends above the main wing, which can increase the negative pressure on the surface of the inner wing segment and increase the lift of the main wing. The leading edge of the wing also generates a large range of vortex, which is attached to the upper surface of the wing, and the vortex generated by the wing tip extends to the rear. As the Angle of attack is pulled up, the vortex intensity gradually increases. At t = 0.455 s, the edge flight-off vortex, leading edge vortex and tip vortex converge in the inner section of the main wing and fall off backward at the trailing edge. As the angle of attack continues to increase, from t = 0.455 s to t = 0.6 s, the excessive angle of attack causes the edge wing shedding vortex to break above the wing, resulting in the drop of lift coefficient and drag coefficient at this stage, and the failure of the elevator surface. In addition, the greater the angle of attack, the greater the degree of rupture. It can be seen from the section vorticity and flow diagram that from the initial moment to 0.455 s, the vorticity of each section gradually increases with the increase of the angle of attack; While from 0.455 s to the end of the pull-up, the vorticity of the wing section is large and changes little with the increase of the angle of attack, but the vorticity of the wing section gradually decreases. In addition, in the initial stage of pull-up, the flow line flows from the edge wing to the leading edge of the wing, and with the increase of the angle of attack, the flow line over the edge wing gradually expands to the inner part of the wing, and the flow separation on the wing surface becomes more intense.
3.4. Common Aileron and Elevator Control of Roll Pitch Channel
Based on the virtual flight simulation platform, this section will simulate the aircraft in the rolling and pull-up mode, assuming that the flight height of the aircraft is 17,000 ft, the flight speed is 600 ft per second, and the rudder deviation is 2˚ for ailerons and 4˚ for elevators through the steering gear. Meanwhile, the northeast earth coordinates, Euler Angle, airspeed, air flow angle and angular velocity of the aircraft are taken as outputs. Analyze the trajectory of the aircraft in 10 seconds and the specific situation of these outputs. The flight path is shown in Figure 13.
As depicted in Figure 14, the simulation of the aircraft’s flight mode involves controlling aileron deflection by 2˚ and elevator control deflection by 4˚. When reaching point a, the vehicle ascended to an altitude of 17,040 ft with a roll angle of 0.03˚, pitch angle of 12.4˚, yaw angle of 0.04˚, and an angle of attack at 8.35˚. Upon reaching state point b, the aircraft climbed to an altitude of 17,100 ft while exhibiting a roll angle of −20˚, pitch angle of −4.09˚, yaw angle of −1.82˚, and an angle of attack at −4.1˚. Continuing to state point c saw the vehicle ascend to a height of 17,150 ft with a roll angle measuring −55˚ along with a pitch angle at
−4.8˚ and yaw angle at 1.85˚; additionally displaying an adverse side slip at −2.1˚. When the aircraft reaches the d state point, it ascends to an altitude of 17,130 ft. The roll angle is −50˚, the pitch angle is 2.5˚, and the yaw angle is −5˚. The angle of attack measures 5˚ while the angle of side slip is −2.5˚. These four state points (a, b, c, and d) indicate that the aircraft is in a stage of small amplitude pulling up. Upon reaching the e state point, located at an east coordinate of −145.3 ft, the aircraft climbs further to an altitude of 17280ft. The roll angle decreases to −34.5˚ while the pitch angle increases to 29˚. Additionally, there is a yaw angle of −27.5˚ with an angle of attack measuring 17˚ and no side slip. When reaching the state point, the aircraft’s east coordinate is −490.3 ft, and it ascended to an altitude of 17,660 ft. The roll angle measures −22.3˚, while the pitch angle stands at 29.5˚. Additionally, the yaw angle is recorded as −27˚, with an angle of attack of 7.8˚.
3.5. Rudder Control of the Yaw Channel
The rudder control of the yaw channel was investigated using the same research methodology as described in Section 3.2. Figure 15 illustrates the variation of yaw angle
with rudder deflection, as well as the change in rudder deflection angle over time. The proportional coefficient (kp) of the PID controller was set to −3, while the integral coefficient (ki) was set to −0.2 and the differential coefficient (kd) to −8. Under this control scheme, at time t = 0.24 s, the yaw angle reached 1.21˚ and at t = 0.99 s, reached 4.94˚. At t = 1.74 s, the yaw angle stabilized at 5.01˚ within an error band of ±5˚ from the desired target value, thus maintaining a steady state thereafter.
4. Conclusion
By using numerical virtual flight simulation method and software platform, the longitudinal maneuvering process of F-16 fighter model is numerically simulated. PID controller is established for multiple Angle of attack state points to realize nonlinear controller design of large Angle of attack pull-up process. In addition, the closed-loop attitude response process is numerically simulated, and the
(a)
(b)
Figure 15. Yaw angle change and yaw rudder deviation change during 0 - 5˚ yaw. (a) Angle of yaw; (b) Rudder deflection angle.
![]()
![]()
Figure 16. Enhanced control of the yaw axis through rudder manipulation.
aerodynamic characteristics and flow field characteristics of the aircraft during high Angle of attack pull-up are studied. Then, by disturbing ailerons and elevators to a certain extent, the 12 state parameters of their output are analyzed, and the nonlinear, unsteady, multi-channel coupling virtual flight of the aircraft after the rudder surface disturbance is realized. In the last section, the rudder’s control of the yaw channel is also investigated. The results show that the control parameter significantly improves the rudder’s response speed and control effectiveness to the yaw channel, which lays a solid foundation for the subsequent research of the rudder’s flight attitude in other maneuvers.