Applsci 13 03518 v2

applied
sciences
Article
Short Landing for Flying-Wing Unmanned Aircraft with
Thrust Vector
Huitao Lyu 1,2 , Yong Yin 1,2 , Zheng Gong 1,2 , Yongliang Chen 1,2, *, Yalei Bai 1,2 and Xueqiang Liu 1,2
1 College of Aerospace Engineering, Nanjing University of Aeronautics and Astronautics,

Nanjing 210016, China; lht@nuaa.edu.cn (H.L.)
2 Key Laboratory of Unsteady Aerodynamics and Flow Control, Ministry of Industry and Information
Technology, College of Aerospace Engineering, Nanjing University of Aeronautics and Astronautics,
Nanjing 210016, China
* Correspondence: chenyl79@nuaa.edu.cn
Abstract: The task of achieving a safe and short landing for a flying-wing unmanned aircraft with
a three-bearing-swivel thrust vector is highly challenging. The process is further complicated due
to the need to switch between multiple control modes, while also ensuring the protection of the
flight boundaries from environmental disturbances and model uncertainties to ensure flight safety.
To address this challenge, this paper proposes a short-landing strategy that employs mixed control
using lift fans, thrust vectors, and aerodynamic control surfaces. The extended state observer (ESO) is
integrated into the inner angular rate control and outer sink rate control to account for environmental
disturbances and model uncertainties. To ensure flight safety, the attainable linear and angular
acceleration is calculated through a trim analysis to determine the command value of velocity and
angle of attack during a short landing. Additionally, a flight boundary protection method is employed
which includes an additional command value of the angle of attack, resulting in a higher probability
of a successful landing. This paper provides a detailed description of the short-landing strategy,
including the control objectives for each phase. Finally, a Monte Carlo simulation is conducted to
evaluate the effectiveness and robustness of the short-landing strategy, and the landing accuracy is
assessed using the circular error probability metric.
Keywords: flying-wing unmanned aircraft; ESO; thrust vector control; short landing; flight
boundary protection
Citation: Lyu, H.; Yin, Y.; Gong, Z.;
Chen, Y.; Bai, Y.; Liu, X. Short
Landing for Flying-Wing Unmanned
Aircraft with Thrust Vector. Appl. Sci.
2023, 13, 3518. https://doi.org/ 1. Introduction
10.3390/app13063518 The rapid development of unmanned aerial vehicles (UAVs) and thrust vector control
Academic Editor: Rosario Pecora
technologies has opened up new possibilities for automatic landing operations, including
vertical landing, shipborne rolling vertical landing (SRVL), and short-landing techniques.
Received: 21 February 2023 Among these, the short landing has gained popularity, as it reduces runway reliance by
Revised: 3 March 2023 lowering the touch-down speed and enhances the aircraft’s bring-back capability. However,
Accepted: 6 March 2023 the application of thrust vector control to short landings poses significant challenges.
Published: 9 March 2023
First, the landing process requires careful consideration of the trim problem that
arises from balancing the aerodynamic and direct forces generated by the thrust vector.
Second, a proper control allocation strategy must be developed to account for the varying
Copyright: © 2023 by the authors.
dynamic responses of the actuators involved. Third, a unified control framework that can
Licensee MDPI, Basel, Switzerland. adapt to control-mode switching is crucial for simplifying the controller design. Fourth,
This article is an open access article the control objects for each landing phase must be defined to ensure smooth operation.
distributed under the terms and Finally, the control system must be robust to model uncertainties during the short-landing
conditions of the Creative Commons process. Addressing these challenges is critical to ensure the safety of short landings with
Attribution (CC BY) license (https:// thrust vector control. By fully considering the above five points, a well-designed control
creativecommons.org/licenses/by/ system can facilitate the successful execution of short-landing techniques, bringing new
4.0/). possibilities and benefits to UAV operations.
Appl. Sci. 2023, 13, 3518. https://doi.org/10.3390/app13063518 https://www.mdpi.com/journal/applsci

Appl. Sci. 2023, 13, 3518 2 of 34
Thrust vector control has been extensively researched, with several notable advance-
ments. Integrated thrust-vector (ITV) control technology has been widely used in post-stall
envelope expansion by implementing the engine nozzle [1]. In history, there have been
attempts to apply thrust vectoring for the takeoff and landing of the Harrier aircraft [2,3].
The Lianistry of Defence has led the development of shipborne rolling vertical landing
(SRVL) technology for F-35B carrier aircraft landing, and thrust vector control has been
applied to improve the performance of a damaged fighter [4]. The F-35B is well-known
for its short takeoff and vertical landing (STVL) capabilities [5], which are supported by
thrust vectoring technology. Additionally, the UK Ministry of Defence [6,7] has improved
the payload performance of weapons—as this was previously a weak point for vertical
landings—through the use of thrust vectoring. On 13 October 2018, British pilot Peter
Wilson successfully completed the first SRVL flight [8], landing the F-35B 755 feet above
HMS Queen Elizabeth’s ski-jump deck, and then taxied 175 feet before using the aircraft’s
brakes to bring it to a standstill from 40 knots. Compared to vertical landing, SRVL lands at
a slow forward speed and uses wing lift to supplement the propulsion system. Reference [9]
developed the optimal trajectory transition controller for a thrust-vectored V/STOL aircraft.
In Reference [10], the problem of the short takeoff curve was solved with the optimal tilt
angle optimization. Different from the straight-line approach [11,12] taken by the F-35B in
SRVL, Nicolas et al. [13] used a curved approach with a fixed-wing UAV and evaluated the
predefined spline path, using flight tests. Furthermore, fast transition maneuvers along
curved trajectories can reduce landing times [14].
Flying-wing unmanned aircraft, with their blended wing body layout [15], have a
high lift-to-drag ratio but suffer from bad static heading stability [16] and control coupling,
posing a significant threat to flight safety. To address this, Wang et al. [17] used the
eigenstructure assignment (EA) method to realize three-axis decoupled control, which
considers flying quality specifications during the control law design process. In addition,
nonlinear dynamic inversion control (NDI) is utilized to deal with nonlinear coupling
dynamics, as X-35B did [18]. However, these model-dependent design approach produce
an unsatisfying result in the presence of model uncertainties. Therefore, robust control [19],
incremental nonlinear dynamic inversion (INDI) control [20,21], and adaptive control [22]
have been applied to improve the robustness of modeling uncertainties. Among these,
active disturbance rejection control (ADRC) [23,24] uses disturbance-observer-based control,
treating model uncertainties and environmental disturbances as lumped disturbances. An
extended state observer (ESO) is designed to estimate and compensate for the lumped
disturbance, showing excellent robustness.
The present study proposes a short-landing strategy utilizing thrust vector control,
which features a continuous curved landing approach. To begin with, a detailed aircraft
dynamics model is developed for landing simulation and performance analysis. To enhance
the safety of the landing process, the attainable area of the linear and angular acceleration
is determined based on the trim calculation to determine the reference values for the angle
of attack and velocity during the short landing. Subsequently, a three-axis decoupled
control law is designed using an extended state observer (ESO). The ESO is used to estimate
and compensate for lumped disturbances, thus ensuring robustness against modeling
uncertainties. To address the switching between different control modes during the landing
process, the angular rate is selected as the control object of the inner loop, and the angular
acceleration commands are allocated based on the different dynamics of the actuators. Fur-
thermore, an overall control structure with multimode control is proposed which clarifies
the control objectives for each phase of the landing process. Additionally, a flight-regime
protection method is introduced to enhance landing safety by incorporating an additional
angle of attack-command compensation. Finally, a nominal simulation and a Monte Carlo
simulation are performed to assess the effectiveness and robustness of the proposed control
strategy. The former involves a simulation without model parameters’ perturbation, while
the latter assesses the robustness of the control strategy against modeling uncertainties.
Appl. Sci. 2023, 13, x FOR PEER REVIEW 3 of 35
the proposed control strategy. The former involves a simulation without model parame-
Appl. Sci. 2023, 13, 3518 3 of 34
ters’ perturbation, while the latter assesses the robustness of the control strategy against
modeling uncertainties.
This paper is organized as follows: Section 2 provides a detailed description of the
This paper
configuration of is
theorganized
flying-wingas follows:
UAV, and Section 2 provides
a dynamic modelaisdetailed
developeddescription of the
to analyze
configuration
aircraft’s of the
behavior flying-wing
during UAV, and a process.
the short-landing dynamicSection
model 3ispresents
developed thetotrim
analyze the
analysis
aircraft’s
results behavior
and definesduring the short-landing
the equilibrium points for process. Sectionangle
the reference 3 presents the and
of attack trim velocity,
analysis
resultsare
which andcritical
defines the equilibrium
parameters pointssafe
for ensuring for and
the reference anglelanding.
efficient short of attack and velocity,
Section 4 com-
which are critical parameters for ensuring safe and efficient short
prehensively elaborates on the design of the control system, including the various landing. Section
control4
comprehensively
loops, allocation ofelaborates on the design
control commands, andof thethe control
overall system,
control including
structure. the various
In Section 5, a
controlshort-landing
novel loops, allocation of control
strategy commands,
incorporating and the
lateral andoverall control
vertical structure.
guidance In Section
is proposed. The5,
a novel short-landing
effectiveness strategy
of the proposed incorporating
strategy is evaluatedlateral and vertical
through numericalguidance is proposed.
simulations, which
Thepresented
are effectiveness of the proposed
in Section 6. Finally, strategy
in Section is 7,
evaluated
the mainthrough
findingsnumerical
of the papersimulations,
are sum-
which are presented in Section
marized, and conclusions are drawn. 6. Finally, in Section 7, the main findings of the paper are
summarized, and conclusions are drawn.
2. Aircraft Dynamics Model
2. Aircraft Dynamics Model
2.1.
2.1. Configuration
Configuration Description
Description
Figure
Figure 1 depicts the
1 depicts the flying-wing
flying-wing UAV UAV thatthat was
was utilized
utilized inin this
this study.
study. The
The UAV UAV is is
equipped with two lift fans and a thrust vector mechanism for
equipped with two lift fans and a thrust vector mechanism for pitch and yaw control. pitch and yaw control. The
lift
Thefans, which
lift fans, are powered
which are powered by Li-Po batteries,
by Li-Po are positioned
batteries, at a fixed
are positioned angleangle
at a fixed at theat front
the
of theof
front fuselage. The thrust
the fuselage. vector vector
The thrust mechanism, consisting
mechanism, of a turbojet
consisting engine and
of a turbojet engine a 3BSD
and
nozzle,
a 3BSD is similaristo
nozzle, that used
similar in the
to that usedF-35B [5].F-35B
in the The thrust-to-weight ratio of theratio
[5]. The thrust-to-weight UAV, of de-
the
spite having two lift fans, is approximately 0.69. Although the
UAV, despite having two lift fans, is approximately 0.69. Although the UAV is unable to UAV is unable to perform
vertical
performlanding
verticalas the F-35B
landing does
as the F-35B[25],does
it can accomplish
[25], short takeoff
it can accomplish shortand shortand
takeoff landing
short
(STOSL) by using its lift fans and vectored thrust. The nozzle pitch
landing (STOSL) by using its lift fans and vectored thrust. The nozzle pitch angle ranges angle ranges from 0 to
90 °, while
from 0 to 90 ◦
the, yaw
while angle
the yawranges from
angle −15 tofrom
ranges 15°. The
−15use to 15 ◦
of vectored
. The usethrust is restricted
of vectored thrust to is
the takeoffto
restricted and
thelanding
takeoffphases, thus resulting
and landing phases, in an resulting
thus improvedintakeoff and landing
an improved takeoffperfor-
and
mance.
landingDuring cruise flight,
performance. Duringthe nozzle
cruise remains
flight, level, remains
the nozzle and the level,
rear engine
and thedoorrearisengine
shut,
door is shut, allowing the UAV to function as
allowing the UAV to function as a standard flying-wing aircraft. a standard flying-wing aircraft.
Figure 1. General diagram of the configuration of the subject flyingwing UAV.

Figure 1. General diagram of the configuration of the subject flyingwing UAV.
Figure
Figure 22 displays
displays the
the aerodynamic
aerodynamic control
control surfaces
surfaces utilized
utilized in
in this
this study;
study; they
they consist
consist
of four elevons for controlling pitch and roll, segmented into elevators and
of four elevons for controlling pitch and roll, segmented into elevators and ailerons.ailerons. In
In
addition,
addition, two
two spoiler-slot
spoiler-slot deflectors
deflectors (SSDs)
(SSDs) are
are employed
employed toto govern
govern yaw
yaw motion.
motion. Relevant
Relevant
aircraft
aircraft parameters
parameters are concisely summarized
are concisely summarized in in Table
Table 1.
1.
Appl.
Appl.Sci. 2023,13,
Sci.2023, 13,3518
x FOR PEER REVIEW 44 of
of 34
35
Figure2.2.Top
Figure Topview
viewofofthe
theflying-wing
flying-wingUAV.
UAV.
Table1.1.Physical
Table Physicalparameters
parametersofofthe
thetarget
targetaircraft.
aircraft.
Parameter
Parameter Explanation
Explanation Value Value
m
m Mass Mass 62 kg 62 kg
SS Wing area
Wing area 1.8 m2 1.8 m2
bb WingspanWingspan 3.2 m 3.2 m
Tv max Maximum thrust of the turbojet engine 270 N
T max Maximum thrust of the turbojet engine 270 N
Tf vmax Maximum force per lift fan 75 N
TIf max Maximum force per lift
Moment of inertia around body x-axis fan 72.1 kg·m275 N
xx
IIyyxx Moment
Moment of inertia
of inertia around around body x-axis 63.7 kg·72.1
body y-axis m2 kg·m2
IIzzyy Moment
Moment of inertia
of inertia aroundaround ·m2 kg·m2
body y-axis 112.3 kg63.7
body z-axis
Izz Moment of inertia around body z-axis 112.3 kg·m2
2.2. Lift Fan Model
2.2. Lift Fan Model
During short-landing operations, the lift fans play a crucial role in controlling the
UAV’sDuring
sink rate short-landing operations,
and counteracting the liftpitch
the negative fans moment
play a crucial causedrole in controlling
by the 3BSD nozzle’sthe
UAV’s sink
pitching. rate and counteracting
To maintain the negative
the UAV’s attitude, pitch moment
a combination of lift caused by the
fan control and3BSD noz-
elevator
zle’s pitching.
input is utilized.To maintain the UAV’s attitude, a combination of lift fan control and eleva-
tor input is utilized.
The lift fan system comprises electric motors, as detailed in Table 2. To obtain a high-
fidelity model,
The lift fanthe force comprises
system and torqueelectric
produced by the
motors, as lift
detailedfans are measured
in Table 2. To during
obtain abench
high-
testing.
fidelityThese
model, measurements
the force and are denoted
torque as T f by
produced andthe Q flift
, respectively, and are during
fans are measured calculated as
bench
outlined in [26]:measurements are denoted
testing. These ( as T f and Q f
, respectively, and are calculated
as outlined in [26]: T f = ρCt ω f 2 D4f
(1)
Q f = ρCq ω f 2 D5f
 Tf = ρ Ct ω f D f
2 4
where ω f is the rotation speed of the culvert  fan, ρ is2 the air density, and Ct and Cq are the (1)
Q f = ρ C qω f The
lift fan’s force and torque coefficient, respectively.
D 5f
least-squares estimation method
is employed to estimate the thrust and torque coefficients, using the bench-test data, as
where ωf is the rotation speed of the culvert fan, ρ is the air density, and Ct and Cq are
illustrated in Figure 3a,b. These coefficients are subsequently utilized in the controller
the lift fan’s force and torque coefficient, respectively. The least-squares estimation
design. While the static bench-test data provide valuable information, they are insufficient
method is employed to estimate the thrust and torque coefficients, using the bench-test
to fully characterize the lift fan’s dynamic response. Therefore, the dynamic torque of the
data, as illustrated in Figure 3a,b. These coefficients are subsequently utilized in the con-
lift fan is computed by incorporating instantaneous electric power:
troller design. While the static bench-test data provide valuable information, they are
Appl. Sci. 2023, 13, 3518 5 of 34

UI
Qf p = (2)
ωf
where U and I are the current and voltage, respectively, and Q f p is the calculated torque.
insufficient to fully characterize the lift fan’s dynamic response. Therefore, the dynamic
torque
Table 2. of the liftparameters
Physical fan is computed
of the liftby incorporating instantaneous electric power:
fan.
UI
Parameter Q =
Explanation Value
(2)
mf Mass of blade culvertf fan
fp
ω
0.092 kg
Df Diameter of fan 0.15 m
where U and I are the current and voltage, respectively, and Qfp is the calculated
Irf Inertia moment of blade culvert fan 0.0004 kg·m2
torque. N Pole pairs of the electric motor 6
(a) (b)
(c) (d)
Figure 3.
Figure 3. Bench
Bench test
test results
results for
for lift
lift fans.
fans. (a)
(a) Lift
Lift fan
fan static
static force
force test
test results
results and
and fitting
fitting curve;
curve; the
the
estimated results are Ct = 4.2278. (b) Lift fan static torque test results and fitting curve; the estimated
estimated results are Ct = 4.2278. (b) Lift fan static torque test results and fitting curve; the esti-
results are C = 0.2303. (c) Recorded dynamic electric power change with rotation speed. (d) Calcu-
mated resultsq are Cq = 0.2303. (c) Recorded dynamic electric power change with rotation speed.
lated results of the dynamic lift-fan torque.
(d) Calculated results of the dynamic lift-fan torque.
Table 2. Physical parameters of the lift fan.
As part of the dynamic modeling of the lift fan system, the derivative of the fan’s
rotationParameter
speed can be represented by the following
Explanation model: Value
mf Mass of blade culvert fan 0.092 kg

. min max K f ω f c (t − τd ) − ω f (t) Ir f , − Q f , Q f p + Q f
D
ω ff (t) = Diameter of fan 0.15 m (3)
Irf Ir f
Inertia moment of blade culvert fan 0.0004 kg·m2
where K f isN the control bandwidthPole pairs
of the of the electric
electronic motor (ECS, K f = 306rad/s),
speed controller
τd is the command delay for the rotation speed (τd = 10 ms), and the rate of rotation speed
As part of the dynamic modeling of the lift fan system, the derivative of the fan’s
rotation speed can be represented by the following model:
ω f ( t ) =
( ( ( ) )
min max K f ω fc ( t − τ d ) − ω f ( t ) I rf , −Q f , Q fp + Q f ) (3)
I rf
Appl. Sci. 2023, 13, 3518 6 of 34
is limited by Q f p and Q f . Finally, the moment of the lift fan, M f , is expressed in the body
axis as follows:
Mf = Ir f ωf × Ω + Qf + (CF − CG + BIAS) × Tf (4)
In Equation (4), we have the following:
h iT
ωf = 0 0 ω f
h iT
Ω= p q r
h iT
Tf = 0 0 − T f (5)
h .
iT
Qf = 0 0 Ir f ω f + Q f
h iT
BIAS = bias cos Φ f sin Φ f 0
where bias is the eccentric size of the lift fan (bias = 0.005 m); CF and CG are the positions
of the lift fan and the center of gravity on the body axis, respectively; Φ f is the rotation
angle of the fan, ranging from 0 to 2π radians; and p, q, and r are the UAV’s angular velocity.
2.3. Thrust Vector Model

Table 3 displays the second-order dynamic system model used to represent the thrust
response to throttle input, accounting for position and rate limitations. It should be noted
that positive rotation of the 3BSD nozzle generates negative moments and that the nozzle
is aligned with the x-axis of the body frame.
Table 3. Effector characteristics of the UAV.
Damping Natural Position

Effector Explanation Rate Limit
Ratio Frequency Limit
δe Elevator angle 0.8 30 rad/s [−30◦ , +30◦ ] ±110 ◦ /s
δa Aileron angle 0.8 30 rad/s [−30◦ , +30] ±110 ◦ /s
δLSSD Left SSD deflection angle 0.8 30 rad/s [0◦ , 60◦ ] ±110 ◦ /s
δRSSD Right SSD deflection angle 0.8 30 rad/s [0◦ , 60◦ ] ±110 ◦ /s
δvq 3BSD nozzle pitch angle 1 10 rad/s [0◦ , 90◦ ] ±30 ◦ /s
δvr 3BSD nozzle yaw angle 1 10 rad/s [−15◦ , +15◦ ] ±30 ◦ /s
δvT Turbojet engine throttle 1 1 rad/s [0, 1] 0.3
The thrust vector changes with the rotation of the nozzle. To avoid a singularity at
90 degrees of pitch, the nozzle is first rotated around its pitch axis by an angle, δvr , and
then the new frame is rotated around its yaw axis by an angle, δvq . The rotation matrix
from the nozzle frame to the body frame is MbT .
To prevent a singularity from occurring at 90 degrees of pitch, the thrust vector is
adjusted through the rotation of the nozzle. Specifically, the nozzle is initially rotated
around its pitch axis by an angle of δvq . Following this, the new frame is then rotated
around its yaw axis by an angle of δvr to arrive at the desired orientation. The corresponding
rotation matrix from the nozzle frame to the body frame is denoted as MbT , and it can be
expressed as follows:
  
cos δvr − sin δvr 0 cos δvq 0 sin δvq
MbT =  sin δvr cos δvr 0  0 1 0  (6)
0 0 1 − sin δvq 0 cos δvq
The resulting expression for the thrust vector in the body frame is as follows:
 

Tvx
 
−δvT Tvmax
 δvT Tvmax cos δvq cos δvr
δ T
 Tvy  = MbT  0 =  vT vmax cos δvq sin δvr  (7)

Tvz 0 −δvT Tvmax sin δvq
The resulting expression for the thrust vector in the body frame is as follows:
Tvx  −δ vT Tv max  δ vT Tv max cos δ vq cos δ vr 

     
Tvy  = M bT  0  =  δ vT Tv max cos δ vq sin δ vr  (7)
T   0   −δ T sin δ vq 

Appl. Sci. 2023, 13, 3518  vz   vT v max
7 of 34
To analyze the moments generated by the thrust vector, it is necessary to take the
cross product of the thrust vector and its corresponding lever arm. This yields the follow-
To analyzewhich
ing equation, the moments
represents generated by the generated
the moments thrust vector, by it
theis thrust
necessaryvector:to take the cross
product of the thrust vector and its corresponding lever arm. This yields the following
equation, which represents M vx  the 0 moments
0  Tvx   by the thrust
0 generated 0 vector: 
      
  Mvy  = 0 0 − Rv  Tvy  = − Rvδ vT Tv max sin δ vq   (8)
Mvx  M 0 00 R 0 0 Tvx T   R δ T 0 

δ δ
 Mvy  = vz0 0 − Rv  Tvy  vz   v −vTRvv δmax Tvmaxvq sin δvq  
v cos sin vr 
= vT  (8)
Mvz
where Rv represents Rv
the0distance 0
from Tvznozzle R
the vT Tcenter
tov δthe vmax cos vq sin δvr
of δgravity along the body’s
x-axis, as shown in Figure 4. To balance the pitch moment generated by the thrust vector,
where Rv represents
the force that the liftthe distance
fans shouldfrom providethe nozzle to the as
is calculated center of gravity along the body’s
follows:
x-axis, as shown in Figure 4. To balance the pitch moment generated by the thrust vector,
the force that the lift fans should provide is Rvδcalculated
T sin δasvq follows:
Tf =
vT v max
(9)
Rv δvT Tvmax nR fsin δvq
Tf = (9)
nR
where n is the number of lift fans (n = 2), and Rf f is the distance between the lift fan center
where
and CG n is the number
along the bodyofx-axis.
lift fans (n = 2),(9),
Equation Rf is the
andwhich distancethe
represents between the lift fan
force required center
to balance
and
the CG
pitchalong the body
moment, x-axis.
is used Equation (9),
to determine thewhich
baseline represents
rotation thespeedforceforrequired to balance
pitch attitude con-
the pitch
trol, moment,
as shown is used
in the to determine
following equation: the baseline rotation speed for pitch attitude control,
as shown in the following equation:
s Rvδ vT Tv max sin δ vq
ω f =Rv δvT Tvmax sin4 δvq (10)
ωf = nR f ρ C t D4 f (10)
nR f ρC D
t f
Figure 4. Longitudinal force analysis.

Figure 4. Longitudinal force analysis.
2.4. Aerodynamic Model
Computational
2.4. Aerodynamic fluid dynamics (CFD) methods are used to determine the UAV’s
Model
aerodynamic
Computational fluidHowever,
properties. dynamicsfor the sake
(CFD) of computational
methods efficiency,the
are used to determine theUAV’s
effectsaer-
of
the lift fans and thrust vector are not considered. Figure 5 displays the pressure
odynamic properties. However, for the sake of computational efficiency, the effects of thecontour of
the aircraft, and three different grids are
lift fans and thrust vector are not considered.employed to check the mesh sensitivity of
Figure 5 displays the pressure contour the CFD
of the
calculation:
aircraft, anda three
fine mesh
different 3.1 ×are
with grids 106employed
nodes, a medium
to check mesh
the 2.2 × 106 nodes,
with sensitivity
mesh of the and
CFD
a coarse mesh with 8 × 105 nodes. In6Figure 6a, the lift coefficient, CL , increases 6 with an
calculation: a fine mesh with 3.1 × 10 nodes, a medium mesh with 2.2 × 10 nodes, and
increase in the angle of attack, but it does not reach the stall region within 12 degrees. The
a coarse mesh with 8 × 105 nodes. In Figure 6a, the lift coefficient, CL , increases with an
drag coefficient, CD , on the other hand, rises monotonically as the angle of attack increases,
increase in the angle of attack, but it does not reach the stall region within 12 degrees. The
as depicted in Figure 6b. It is worth noting that the drag coefficient obtained using the
drag coefficient, CD , on the other hand, rises monotonically as the angle of attack in-
medium mesh is quite close to that obtained by the fine mesh. Therefore, unless otherwise
creases, as depicted in Figure 6b. It is worth noting that the drag coefficient obtained using
stated in subsequent discussions, most of the results reported below are obtained on the
the medium mesh is quite close to that obtained by the fine mesh. Therefore, unless oth-
medium mesh.
erwise stated in subsequent discussions, most of the results reported below are obtained
The pitch moment coefficient, Cm , in Figure 6c indicates that the UAV gains longitu-
on the medium mesh.
dinal static stability before 10 degrees, so the maximum angle of attack for safe flight is
set to 10 degrees. However, the UAV exhibits yaw static instability with a negative static
directional-stability derivative Cnβ < 0 [27], which makes yaw control more challenging.
This is shown in Figure 6d.
Appl.
Appl.Sci. 2023,13,
Sci.2023, 13,3518
x FOR PEER REVIEW 88 of
of 34
35
Figure 5. Pressure contour of the aircraft.
The pitch moment coefficient, Cm , in Figure 6c indicates that the UAV gains
dinal static stability before 10 degrees, so the maximum angle of attack for safe flig
to 10 degrees. However, the UAV exhibits yaw static instability with a negative s
rectional-stability derivative Cnβ < 0 [27], which makes yaw control more chall
Figure5.5.Pressure
Figure Pressurecontour
contourof
ofthe
theaircraft.
aircraft.
The pitch moment coefficient, Cm , in Figure 6c indicates that the UAV gains longitu-
dinal static stability before 10 degrees, so the maximum angle of attack for safe flight is set
to 10 degrees. However, the UAV exhibits yaw static instability with a negative static di-
rectional-stability derivative Cnβ < 0 [27], which makes yaw control more challenging.
(a) (b)
(a) (b)
(c) (d)
Figureaerodynamic
Figure 6. Longitudinal 6. Longitudinal aerodynamic
coefficients: (a) liftcoefficients: (a) lift coefficient
coefficient (sideslip angle β =(sideslip angle β = 0),
0), (b) drag
coefficient (sideslip angle β = 0), (c) pitch-moment coefficient, and (d) directional static stability static sta
coefficient (sideslip angle β = 0), (c) pitch-moment coefficient, and (d) directional
rivative (per radian).
derivative (per radian).
The SSD is used to The SSD is

control theused
yawtomotion
control ofthe
the yaw
UAV.motion
However,of the
theUAV.
controlHowever,
derivative, the contro
CnδSDD(c) ative, C ,
, exhibits nonlinear behavior
nδ exhibits nonlinear behavior
at deflection angles below 10 degrees due to interactions degrees
at deflection angles below 10
SDD (d)
interactions between the SSD and
between the SSD and the wing’s boundary layer [28], as shown the wing’sinboundary layernonlinearity
Figure 7. This [28], as shown in F
Figure 6. Longitudinal
This aerodynamic
nonlinearity coefficients:
makes yaw (a) lift
control coefficient
more (sideslip
challenging.
makes yaw control more challenging. To address this, a bias angle [29] of 10 degrees angle
To β = 0), (b)this,
address drag
is a bias an
coefficient (sideslip angle
of 10 β =
degrees 0), (c) pitch-moment coefficient, and (d) directional static stability de-
chosen for the SSD during yawiscontrol.
chosen for the SSD during yaw control.
rivative (per radian).
The SSD is used to control the yaw motion of the UAV. However, the control deriv-
ative, CnδSDD , exhibits nonlinear behavior at deflection angles below 10 degrees due to
interactions between the SSD and the wing’s boundary layer [28], as shown in Figure 7.
This nonlinearity makes yaw control more challenging. To address this, a bias angle [29]
of 10 degrees is chosen for the SSD during yaw control.
Appl. Sci. 2023, 13, x FOR PEER REVIEW
Appl. Sci. 2023, 13, 3518 9 of 34
Figure
Figure 7.7.SSD
SSD yaw
yaw control
control moment
moment derivative
derivative (perthe
(per radian); radian); the
gray area is gray area is the
the continuous continuous
error bar e
for various
for variousangles of attack
angles (angles
of attack of attack
(angles = −6~12
of αattack = 6~12 deg).
α deg).
Equations of motion with six degrees of freedom are developed for the UAV based on
Equations
the previous of motion
analysis. with six
These equations aredegrees
then usedofforfreedom are developed
the performance for the UAV
analysis, control
on the design,
system previous
andanalysis.
simulationThese
of shortequations are then
landings. Table used for
3 summarizes thethe performance
simulation inputs, analys
and the
trol outputs
system are the and
design, UAV’s states duringofsimulation,
simulation includingTable
short landings. position, attitude angles,the sim
3 summarizes
and angular rates.
inputs, and the outputs are the UAV’s states during simulation, including positio
tude angles, andPerformance
3. Short-Landing angular rates.
Analysis
3.1. Trim for Short Landing
3. Short-Landing
To land the UAVPerformance Analysis
safely, the sink rate, angle of attack, and velocity must be maintained
for a steady approach. However, determining the baseline values for these parameters is
3.1. Trim for Short Landing
challenging. To address this, the target UAV is trimmed at various angles of attack, sink
rates,To
andland the UAV
velocities safely,
at a height of 50the sink rate,
m under angle
control input of attack, and
restrictions. The velocity
trim speedmust be
ranges from 17 to 50 m/s, the angle of attack from 3 to 12 degrees,
tained for a steady approach. However, determining the baseline values and the target sink rate
for these
from 0 to − 4 m/s. The lift fan control is split into two parts: balancing
eters is challenging. To address this, the target UAV is trimmed at various the moment fromangles of
δvq and compensating for δe , as shown in Figure 8.
sink The
rates, andinvelocities
results Figure 9 showat a that
height of 50 m area
the reachable under control
decreases as input
the anglerestrictions.
of at- T
speed ranges and
tack increases, fromthe17velocity
to 50 m/s, the angle
boundaries moveoftoattack from
the left. As 3 thetoangle
12 degrees,
of attack and the
sink rate from
increases, 0 to boundary
the speed −4 m/s. The lift fanallowing
decreases, control for
is split intoflight
a lower twospeed
parts:for
balancing
the ap- the m
proach. However, a lower speed comes with a reduced
from δvq and compensating for δe , as shown in Figure 8. maneuvering area, which increases
the risk of the approach phase. Therefore, a balance must be made between the land-
ing speed and flight safety, and the proper angle of attack and flight velocity must be
selected during the short landing.
Appl. Sci.
Appl. Sci. 2023,
2023, 13,
13, 3518
x FOR PEER REVIEW 10
10 of 35
of 34
Figure 8.
Figure 8. Trim
Trim procedure
procedure for
forshort-landing
short-landingperformance
performanceanalysis.
analysis.
The results in Figure 9 show that the reachable area decreases as the angle of attack
increases, and the velocity boundaries move to the left. As the angle of attack increases,
the speed boundary decreases, allowing for a lower flight speed for the approach. How-
ever, a lower speed comes with a reduced maneuvering area, which increases the risk of
the approach phase. Therefore, a balance must be made between the landing speed and
flight safety, and the proper angle of attack and flight velocity must be selected during the
short landing.
Figure 9. Reachable area for different angles of attack (angle of attack, α = 3~12 deg).
Figure 9. Reachable area for different angles of attack (angle of attack, α = 3~12 deg).
3.2. Attainable Linear and Pitch Angular Acceleration

To ensure flight safety during short landings, it is important to select an appropriate
Appl. Sci. 2023, 13, 3518 11 of 34
3.2. Attainable Linear and Pitch Angular Acceleration

To ensure flight safety during short landings, it is important to select an appropriate
angle of attack and velocity for the aircraft. The conventional flying quality requirements
based on maximum angular acceleration and attitude change in 1 s under an abrupt
(0.3-s) ramp-control input [30] are not suitable for the flying-wing UAV, due to its unique
size and inertial properties. Instead, a methodology can be used to identify the safest angle
of attack for a given velocity by evaluating the maximum linear and angular acceleration
generated by control effectors. This methodology involves selecting a series of equilibrium
points [31] for the aircraft during a short landing. The longitudinal dynamic model of the
UAV can be derived from Figure 4.
.
V gx = − m1 D cos γ + L sin γ − T cos δvq + θ + Tv RRv sin δvq sin θ

f
. (11)
V gz = − m − D sin γ + L cos γ + T sin δvq + θ − gm + Tv RRv sin δvq cos θ
1

f
Equation (11) describes the longitudinal dynamic model of the UAV during a short
landing, where Vgx and Vgz are the velocity components in the earth coordinate frame;
L and D are the lift and drag forces in the wind coordinate frame; γ and θ are the flight
path angle and pitch attitude angle, respectively; and g is the gravitational acceleration. To
simplify the analysis, Equation (11) is expanded using a Taylor series, resulting in Equation
(12), which considers only the first-order term of the expansion:
. . ∂ f ( x, u) ∂ f ( x, u)
x ≈ x0 + ( x − x0 ) + ( u − u0 ) (12)
∂x x = x0 ∂u x = x0
u = u0 u = u0
| {z } | {z }
F (·) G (·)
where the state vector is x = [Vgx Vgz ]T , and the control input vector is u = [δvq Tv ]T ; x0 and
u0 are the trimmed values for state and control; and F(·)∈R2×2 and G(·)∈R2×2 are vectors
of nonlinear algebraic equations. The control input matrix, G(·), represents the acceleration
response to the control inputs, as shown in the following equation:
 . . 
∂V gx ∂V gx
 ∂δ ∂Tv 
G (·) =  . vq . (13)
∂V ∂V gz

gz
∂δvq ∂Tv
Equation (13) can be expanded using partial derivatives, as demonstrated in

Equation (14).
. .
∂V gx ∂V gx
= − Tmv sin δvq + θ + RRvf cos δvq sin θ
1
Rv
∂δvq ∂Tv = m cos δvq + θ − Rf sin δvq sin θ
. . (14)
∂V gz Tv
Rv ∂V gz 1
Rv
∂δvq = − m cos δvq + θ + R cos δvq cos θ
f ∂Tv = − m sin δvq + θ + R sin δvq cos θ
f
Figure 10 displays the matrix G(·) calculated at a trimmed angle of attack of 4 degrees
and various velocities. The gray area indicates a continuous error bar for the changing sink
rate. In Figure 10a, tangential acceleration shows a negative correlation with the nozzle
pitch angle, indicating that the UAV gains more deceleration capability as velocity decreases
due to the larger trimmed nozzle pitch angle. In contrast, the throttle shows a positive
correlation with tangential acceleration, which gradually weakens as the velocity decreases
in Figure 10b. The nozzle pitch angle and throttle both exhibit a negative correlation
with the normal acceleration in Figure 10c,d. The sink rate significantly affects the control
efficiency for normal acceleration when using the nozzle pitch angle. The throttle shows
better control for normal acceleration than the nozzle pitch angle. Overall, for better control
effectiveness, it is preferable to use the throttle to control Vgz and the nozzle pitch angle to
control Vgx , despite the control coupling phenomenon between the two.
decreases in Figure 10b. The nozzle pitch angle and throttle both exhibit a negative corre-
lation with the normal acceleration in Figure 10c,d. The sink rate significantly affects the
control efficiency for normal acceleration when using the nozzle pitch angle. The throttle
shows better control for normal acceleration than the nozzle pitch angle. Overall, for better
Appl. Sci. 2023, 13, 3518 control effectiveness, it is preferable to use the throttle to control Vgz and the nozzle
12 ofpitch
34
angle to control Vgx , despite the control coupling phenomenon between the two.
(a) (b)
(c) (d)
Figure
Figure 10.10.Control
Controlinput
inputmatrix
matrixfor
for acceleration
acceleration (angle
(angle of
of attack
attackαα==4 4deg);
deg);the gray
the grayarea is the
area con-
is the
tinuous error bar for various sink rates (0~−4m/s): (a) G11, (b) G12, (c) G21, and (d) G22.
continuous error bar for various sink rates (0~−4m/s): (a) G11 , (b) G12 , (c) G21 , and (d) G22 .
Equation (13) allows the acceleration control matrix to be calculated from trimmed
results at different velocities and angles of attack. The attainable acceleration set (AAS)
area can be evaluated from this matrix. The AAS area decreases with increasing trimmed
velocity and increases with decreasing absolute sink rate, as shown in Figure 11a. In
Figure 11b, the AAS area and attainable pitch angular acceleration (APAA) are compared,
and they move in opposite directions as velocity increases. The equilibrium point, where
the normalized AAS and APAA curves intersect, is shown in Figure 11c for different sink
rates. The equilibrium points are calculated at various angles of attack and sink rates in
Figure 11d, Figure 12a shows the AAS area [32] for a single trimmed point, with the throttle
ranging from −10% to 10% and nozzle pitch angle from −10 to 10 degrees at an angle of
attack of 6 degrees, and their mean velocity is presented in Figure 12b with an error bar.
The equilibrium curve between the velocity and angle of attack provides the appropriate
control command value for flight safety. Therefore, landing at the angle of attack and speed
in Figure 12b maximizes the maneuvering margin during landing.
For flight safety during short landings, Figure 12b indicates an entering speed of
45 m/s in level flight condition with the nozzle pitch angle down approximately 45 degrees,
as illustrated in Figure 13a. The lift fan should be set at 20%, as shown in Figure 13b. For
continuous approaches, the recommended command values are a velocity of 35 m/s and
angle of attack of 6 degrees, accounting for aerodynamic control effectiveness requirements
for disturbance rejection.
Figure 11d, Figure 12a shows the AAS area [32] for a single trimmed point, with the throt-
tle ranging from −10% to 10% and nozzle pitch angle from −10 to 10 degrees at an angle of
attack of 6 degrees, and their mean velocity is presented in Figure 12b with an error bar.
The equilibrium curve between the velocity and angle of attack provides the appropriate
Appl. Sci. 2023, 13, 3518 control command value for flight safety. Therefore, landing at the angle of attack13andof 34
speed in Figure 12b maximizes the maneuvering margin during landing.
(a) (b)
(c) (d)
Figure 11. AAS and AMS calculation for short landing (angle of attack = 6 deg): (a) AAS area results,
Figure 11. AAS and AMS calculation for short landing (angle of attack = 6 deg): (a) AAS area results,
(b) AAS area and AMS results, (c) normalized AAS area and AMS results, and (d) equilibrium points
(b) AAS area and AMS results, (c) normalized AAS area and AMS results, and (d) equilibrium points
for different sink rates.
for different sink rates.
(a) (b)
Figure
Figure 12.12. AAS
AAS area
area and
andequilibrium
equilibriumpoints curve:
points (a) AAS
curve: area (angle
(a) AAS of attack
area (angle of =attack
6 deg,=sink rate =
6 deg,
−1 m/s, and velocity = 35 m/s) and (b) equilibrium points curve.
sink rate = −1 m/s, and velocity = 35 m/s) and (b) equilibrium points curve.
For flight safety during short landings, Figure 12b indicates an entering speed of 45
m/s in level flight condition with the nozzle pitch angle down approximately 45 degrees,
angle of attack of 6 degrees, accounting for aerodynamic control effectiveness require-
ments for disturbance rejection.
For flight safety during short landings, Figure 12b indicates an entering speed of 45
m/s in level flight condition with the nozzle pitch angle down approximately 45 degrees,
Appl. Sci. 2023, 13, 3518 angle of attack of 6 degrees, accounting for aerodynamic control effectiveness require-
14 of 34
ments for disturbance rejection.
(a) (b)
Figure13.
Figure 13.Trim
Trimresults
resultsfor
for short
short landing
landing (angle
(angle ofof attack
attack = =3 3deg):
deg):(a)
(a)throttle
throttleand
andnozzle
nozzleangle
angletrim
trim
results and (b) throttle and lift fan trim results.
results and (b) throttle and lift fan trim results.
4.4.Flight
FlightControl
ControlSystem
SystemDesign
Design
4.1.
4.1.Angular
AngularRateRateand
andAttitude
AttitudeController
Controller
InInananinner
innerangular
angularrate
ratecontroller
controllerloop,
loop,the
therelationship
relationshipbetween
betweenthe
thecontrol
controlinput
input
and
andthethedynamics
dynamicsofofthetheUAV
UAV is is
obtained
obtainedbyby
the following
the following equation:
equation:
.
Ω = I −Ω
1 M
= (·)
−1 1
M(I⋅) − I{−Ω
I −1 − {Ω××IΩ
I Ω}}+
+ vv((·) (15)
⋅)
| 
{z } (15)
f (·) f(⋅)
whereΩΩ==[p,
where q, q,
[ p, r]Tr]represents
T represents the
theroll,
roll,pitch,
pitch,and
andyaw
yawangular
angularrates
ratesininthe
thebody
bodyframe,
frame,
T
respectively;vv(·(·)) =
respectively; = [v[vpp, ,vqv,qv, r]vTr is
] the
is the control
control acceleration
acceleration input;
input; I ismoment
I is the the momentof the of
inertialthe
inertial
matrix;matrix;
and f(·)and
is the f (·)total
is the total disturbance
disturbance added added to the flying-wing
to the flying-wing UAV for UAV foraxis,
each eachwhich axis,
which is composed
is composed of modeled
of modeled couplingdynamics
coupling dynamicsand and unmodeled dynamics, MM
unmodeled dynamics, (·)(·). .The
The
unmodeled dynamics include external environmental disturbance and internal coupling
dynamics. Then the Equation (15) can be reformulated as follows:
( .
x = f (·) + v(·)
. (16)
f (·) = w(t)
To control the angular acceleration around the UAV’s body axis, an ESO-based control
law is designed, where w(t) is the derivative of f (·) that is unknown but bounded. An ESO is
capable of estimating the lumped disturbance caused by external disturbances, unmodeled
dynamics, parameter variations, and complex nonlinear dynamics [33]. Dynamic compen-
sation is then used based on the estimated results to linearize the model. For roll, pitch,
and yaw control, a second-order linear ESO is designed, and the observer is expressed as
follows:
.
z1 = β 01 (y − z1 ) + z2 + v(·)
. (17)
z2 = β 02 (y − z1 )
In each angular rate control loop, y represents the angular rate measurements, while
z1 and z2 are the estimated angular rates and disturbance. The observer gains, β01 and β02 ,
can be determined using bandwidth methods [34].
β 01 = 2ωo,(·) , β 02 = ωo,(·) 2 (18)
The ESO’s observation bandwidth for each angular rate control loop is denoted as
ωo,(·) . The error dynamics of the ESO can be calculated as follows:
Appl. Sci. 2023, 13, 3518 15 of 34
.
e1 = − β 01 e1 + e2
. (19)
e2 = − β 02 e1 − w(t)
where e1 = z1 − x and e2 = z2 − f (·) are the estimated error. The characteristic polynomial
of Equation (19) can be expressed as follows:
λ(s) = s2 + β 01 s + β 02 (20)
With positive observer gains, the roots of the characteristic polynomial in Equation (20)
are all in the left half plane, ensuring the ESO’s bounded-input–bounded-output (BIBO)
stability [35]. A proportional controller is designed, and the estimated disturbance is given
as follows:
v(·) = k (·) (Ωc − Ω) − z2(·) (21)
The proportional controller is designed to control the estimated disturbance in each
control channel. Theproportional gains or control bandwidth of the closed loop are denoted
by k(·) ∈ (k p , kq , kr . The estimated disturbance for each control channel is denoted by z2(·) ,
and Ωc = [ pc , qc , rc ]T is the command angular rate. The stability analysis of the ESO-based
controller and the relationship between the error bounds and ADRC bandwidth are presented
in Reference [34]. The rate of change in attitude angles can be expressed as follows:
.   
φ
.
1 sin φ tan θ cos φ tan θ p
θ = 0 cos φ − sin φ  q  (22)
  
.
ψ 0 sin φ/cos θ cos φ/cos θ r
where φ, θ, and ψ are roll, pitch, and yaw angles, respectively. The desired rotational rates
can be achieved using Equation (22), which gives the following [36]:
−1
pc 1 sin φ tan θ vφ tan θ/cos φ
= − r (23)
qc 0 cos φ vθ − tan φ
where vφ and vθ are the desired dynamic responses for roll and pitch attitude, which are
calculated using the following equation:

vf k f 0 f c − f pc
= (24)
vθ 0 k θ θc − θ qc
The slow outer attitude control loop [37] is governed by an NDI-based control law,
which is a model-based strategy for controlling nonlinear dynamic systems. The input-
output linearization method, based on NDI control theory [38], is utilized due to the precise
measurement of the aircraft’s kinematic parameters, using an inertial navigation system
(INS). This enables the design of a simple linear controller for the resulting system, where
the command attitude angles [φc , θ c ] and proportional control gains [kφ , kθ ] are utilized. To
achieve coordinated flight, it is necessary to regulate the sideslip angle, β, instead of the
yaw angle ψ, particularly for tailless aircraft. The desired yaw angular rate can be obtained
by using the following equation, as proposed in Reference [39]:
g sin φ cos θ
rc = −k β ( β c − β) + (25)
V
where βc is the command sideslip angle, which is set to zero throughout the short-landing
process, and kβ is the control gain. The complete diagram for the attitude and angular rate
control loop is illustrated in Figure 14.
obtained by using the following equation, as proposed in Reference [39]:
g sin φ cos θ
rc = − kβ ( β c − β ) + (25)
V
where βc is the command sideslip angle, which is set to zero throughout the short-landing
Appl. Sci. 2023, 13, 3518 16 of 34
process, and kβ is the control gain. The complete diagram for the attitude and angular rate
control loop is illustrated in Figure 14.
Figure 14. Complete control diagram for inner and outer

outer loop
loop control.
control.
4.2. Control Allocation

4.2. Control Allocation
In the UAV configuration depicted in Figure 1, both the roll and yaw control are
In the UAV configuration depicted in Figure 1, both the roll and yaw control are
achieved by the aileron and SSD. The control allocation process in this configuration can
achieved by the aileron and SSD. The control allocation process in this configuration can
utilize a pseudo-inverse method, which is outlined below:
utilize a pseudo-inverse method, which is outlined below:
−1
00  −1  vv p

δac δ  LδLa
= δ (26)
 = 0 NδSDD   vr
ac p

a
δSDDc (26)
δ SDDc   0 NδSDD   vr 
where δac and δSDDc are the aileron and SSD command values, and the roll and yaw control
where δac andparameters
effectiveness δSDDc are theare
aileron and SSD command values, and the roll and yaw control
as follows:
effectiveness parameters are as follows:1/2ρV 2 SbC
lδa
Lδa = I
Appl. Sci. 2023, 13, x FOR PEER REVIEW ρxxV 22SbC
1 21/2ρV lδ a
17 of(27)
35
NδLSDD
δa
== SbCnδ
SDD
I xx Izz
(27)
where Clδa is the control derivative of the aileron. 1 2 ρ VThe 2
SbC control
nδSDD allocation, Equation (26), does
disturbance that can be estimated by Nthe =
ESO. The control allocation for the pitch channel
δSDDthe SDD and aileron, as it is treated as a disturbance
not consider the coupling effect between I zz from the roll and yaw channels. Due
involves
that can beboth aero surfaces
estimated by theand
ESO.liftThe
fans, which
control differ
allocation for the pitch channel involves both
to thesurfaces
where
aero distinct response
Clδa is and
the control characteristics
lift fans,derivative
which differ of of theaileron.
the
from elevator
the roll The
and and lift channels.
control
yaw fans, bothDue
a high-frequency
allocation, Equation (26),
to the distinct
and low-frequency
does not characteristics
response control
consider the couplingallocation
effectand
of the elevator for the
between pitch
lift fans, the channel
SDD
both are carried out,
and aileron, asand
a high-frequency astreated
it is illustrated
low-frequencyas a
in Figure
control 15.
allocation for the pitch channel are carried out, as illustrated in Figure 15.
Figure 15.
Figure Pitch channel
15. Pitch channel control
control allocation.
allocation.
To allocate the control inputs for the pitch channel, the desired pitch angular accel-
To allocate the control inputs for the pitch channel, the desired pitch angular accel-
eration, vq , is first filtered using a first-order lag transfer function with a time constant of
eration, vq, is first filtered using a first-order lag transfer function with a time constant of
1 s (τ = 1 s). The resulting desired pitch control value is then split into two components: a
1 s (τrr = 1 s). The resulting desired pitch control value is then split into two components:
low-frequency part (LFP), v ; and a high-frequency part (HFP), vqh . The LFP is allocated
a low-frequency part (LFP), ql vql; and a high-frequency part (HFP), vqh. The LFP is allocated
using a daisy-chaining method [32], which involves calculating the elevator command
using a daisy-chaining method [32], which involves calculating the elevator command
value as follows:
value as follows:

δecs = max min vql /Mδe , δemax , δemin (28)
( ( )
δ ecs = max min vql Mδ , δ e max , δ e min
e
) (28)
where Mδe is the control effectiveness parameter of the elevator and can be calculated as
follows:
1 2 ρ V 2ScCmδe
Mδe = (29)
I yy
Appl. Sci. 2023, 13, 3518 17 of 34
where Mδe is the control effectiveness parameter of the elevator and can be calculated as follows:
1/2ρV 2 ScCmδe
Mδe = (29)
Iyy
where Cmδe is the control derivative of the elevator. If the elevator deflection angle reaches
its limits, the lift fan is utilized to compensate for the LFP. The left LFP, uqs , is then obtained
through the lift fan:
vqs = vql − Mδe δecs
(30)
ω f cs = vqs /M f an
where ω fcs is the command value of the lift fan, and Mfan is the control effectiveness
parameter of the lift fan, as follows:
2nρCt ω f D4f R f
M f an = (31)
Iyy
Secondly, the HFP is allocated via parallel allocation. The elevator’s and lift fan’s
residual control powers are denoted as urde and urfan , respectively. When vqh ≥ 0, the output
of the parallel output is given by the following:


 k = v qh / u rde + u r f an



ω f cp = k ω f max − ω f cs (32)


 δecp = k(δemax − δecs )


otherwise, we have the following:



 k = v qh / u rde + u r f an



ω f cp = k ω f min − ω f cs (33)


 δecp = k(δemin − δecs )


where we have the following:
 Mδe (δemax − δecs ), vqh < 0


urde =
 Mδe (δemin − δecs ), vqh ≥ 0
 (34)

 M f an ω f max − ω f cs , vqh ≥ 0
ur f an =
 M f an ω f min − ω f cs , vqh < 0

In Equation (34), ω f max and ω f min represent the maximum and minimum rotation
speeds of the lift fans, respectively. These values are selected to be between −80% and 80%
of the baseline rotation speed command, which is calculated using the throttle command
value in Equation (10). To synchronize the responses of the turbojet engine and lift fans,
a first-order lead-lag controller is applied with time constants of 0.1 s (τ f = 0.1 s) and 1 s
(τt = 1 s) for the engine and lift fans, respectively. The baseline rotation speed command
passes through the lead-lag control unit and is denoted as ω f cb in Figure 15. The output
of the pitch control allocation is a combination of the baseline part, the daisy chaining
allocation part, and the parallel allocation part, as expressed by the following equation:

ω f c = ω f cb + ω f cs + ω f cp
(35)
δec = δecs + δecp
Appl. Sci. 2023, 13, 3518 18 of 34
4.3. Angle of Attack Control

For flying wing aircraft, changes in the angle of attack significantly affect longitudinal
static stability. Poor longitudinal and lateral static stability can pose a serious threat to
flight safety. Therefore, the angle of attack, denoted by α, is chosen as the control object
for the outer loop, which generates a pitch angular rate command for the inner control
loop. Neglecting the lateral and directional coupling dynamics, the linearized model for
the angle of attack can be expressed as follows [17]:
.
α = − Zα α + q + Zδe δe (36)
1/2ρV 2 SCLα
Zα = mV
Appl. Sci. 2023, 13, x FOR PEER REVIEW 19 of(37)
35
1/2ρV 2 SCLδe
Zδe = mV
where CLα is the lift curve slope, and CLδe is the derivative of the lift coefficient to the eleva-
α ( s)
tor. The value of CLδe is insignificant when compared to CLα , resulting in the simplification
1
of the transfer function for the pitch angular=rate to angle of attack: (38)
q ( s ) s + Zα
α(s) 1
= (38)
Consider the closed-loop dynamics q(s)of the s +pitch
Zα angular rate in the system. If it is
assumed to be athe
Consider first-order low-pass
closed-loop filter, then
dynamics thepitch
of the transfer function
angular ratefor
in the
the angle of attack
system. If it is
command
assumed to be a first-order low-pass filter, then the transfer function for the angle ofmodel
to the measured output can be simplified as a second-order nominal attack
with the angle
command of measured
to the attack controller. Thisbe
output can can be expressed
simplified as follows: nominal model with
as a second-order
α ( s)
the angle of attack controller. This can be expressed as follows:
kα kq
= 2 kα kq (39)
α c ( s=
α(s)
) s2s ++kkqss ++ kkα kkq (39)
αc (s) q α q
The
The control diagram of the angle of attack is shownin
control diagram of the angle of attack is shown inFigure
Figure16.
16.
Figure16.
Figure Diagramofofthe
16.Diagram theangle
angleofofattack
attackcontroller.
controller.
The controller gain, k , can be calculated using the desired bandwidth and damping
The controller gain, kαα, can be calculated using the desired bandwidth and damping
ratio of the closed-loop system. In this case, the damping ratio is selected to be 0.8, and
ratio of the closed-loop system. In this case, the damping ratio is selected to be 0.8, and
then kα can be calculated using the following formula:
then kα can be calculated using the following formula:
k α = k q /1.622 (40)
kα = kq 1.6 (40)
4.4. Sink Rate and Velocity Control
According
4.4. Sink Rate and to the analysis
Velocity Control in Section 3, the 3BSD’s nozzle pitch angle is utilized to
control the UAV’s velocity, and its sink rate is controlled by using the turbojet engine
According to the analysis in Section 3, the 3BSD’s nozzle pitch angle is utilized to
throttle. For velocity control, a simple proportional-integral (PI) controller is applied to
control the desired
obtain the UAV’s tangential
velocity, and its sink rate
acceleration, is controlled
as the by using the turbojet engine
following equation:
throttle. For velocity control, a simple
Z t proportional-integral (PI) controller is applied to
.
obtain the desired tangential acceleration,
k as the following
(Vc − V )dt − k vp V = V equation: (41)
vi des
0
kvi  ( Vc − V ) dt − k vp V = V des
t
where we have the following: 0

(41)
k vp = 2ζ T ωT
where we have the following: (42)
k vi = ωT2
kvp = 2ζ T ωT
(42)
kvi = ωT2
where ωT is the closed-loop bandwidth, and ζT is the damping ratio which is set to 0.8
to reduce overshoot. The command value of the nozzle pitch angle is achieved using the
Appl. Sci. 2023, 13, 3518 19 of 34
where ω T is the closed-loop bandwidth, and ζT is the damping ratio which is set to 0.8
to reduce overshoot. The command value of the nozzle pitch angle is achieved using the
Appl. Sci. 2023, 13, x FOR PEER REVIEW 20
following equation:
  .  
−m V des − g sin γ
= maxminperiod.
δvqc deceleration , δvqmax −
δ dynamics ofδthe
, −δvqb  (43)
δvTcThe closed-loop
cos δvqb vqb vqb pitch angle are regarded as a

Tvmax
order low-pass filter, as shown in Figure 17a. The closed-loop transfer function fro
where δvTc is theto throttle command value, δvqb is the baseline nozzle pitch angle given by
γ is as follows:
the setting at the beginning of short landing, and δvq max is the maximum nozzle pitch
angle. The output command of the nozzle pitch angle γ ( s ) is equalkθtokγ Zthe
α sum of δvqc and δvqb .
= 2
For sink rate control, two controllers are used ( s) different
γ c for s + kθ s +short-landing
kθ kγ Zα stages. The
first uses the pitch angle to control the sink rate [22], and the controller is activated in
To reduce
the deceleration period. The overshoot,
closed-loopthe dampingofratio
dynamics the of the angle
pitch second-order system
are regarded asina Equation
is selected
first-order low-pass filter,toasbeshown
0.8. The in controller
Figure 17a. gain,
Thekclosed-loop
γ, can then be calculated
transfer through
function pole placem
from
kγ = kθ ( 1.6 2 Zα )
γc to γ is as follows:
γ(s) k θ k γ Zα
The other sink rate=controller2 uses the throttle to control the sink rate (44)
and is activ
γc ( s ) s + k θ s + k θ k γ Zα
during the approach stage, as depicted in Figure 17b.
(a)
( b)
Figure
Figure 17. Sink rate 17. Sink
control rate control
diagram: diagram:
(a) control sink (a)
ratecontrol sinkusing
obtained rate obtained usingangle
pitch attitude pitchand
attitude angl
(b) control sink rate obtained using
(b) control sink rate obtained using the throttle. the throttle.
For sink
To reduce overshoot, therate control
damping in the
ratio approach
of the stage, system
second-order the vertical dynamics
in Equation (44)of
is the UAV
given
selected to be 0.8. bycontroller
The the following
gain,equation:
kγ , can then be calculated through pole placement:
α δ vT Tv max sin δ vq nρCtω f D f (45)

ρ 2 2 2 4
k = 1 k 2 / V SC
1.6 Z sin
− V gz cos γ = Zα α V cos α − g cos γ +
γ θ Dα
+ +
m m m
The other sink rate controller uses the throttle to control the sink rate and is activated
where stage,
during the approach we have the following:
as depicted in Figure 17b.
For sink rate control in the approach stage, the vertical dynamics of the UAV are given
by the following equation: V gz = −h
Equation (46) can be rewritten as follows:

. 1/2ρV 2 SCD sin α δvT Tvmax sin δvq nρCt ω f 2 D4f
− V gz /cos γ = Zα αV cos α − g cos γ + + + (46)
m nρ Ct D f Tv max sin δ vq
4
m m
Vgz = −Zα α V + g cos γ − ω fcb 2 − δ vTc + fa ( ⋅)
where we have the following:  m    m 
. .
A b
V gz = −h 1 0
(47)
where
Equation (46) canAbe
1 is the known part of the sink rate control, b0 is the control effectiveness param
rewritten as follows:
for the throttle, and fa(·) is the lumped disturbance, which includes the unmodeled
namics and external disturbances. An NDI controller is designed as a baseline sink
controller, as shown in Figure 17b. The output of the baseline controller is formulat
follows:
( )
azE ,NDI = − ksr hc − h − A1
Appl. Sci. 2023, 13, 3518 20 of 34
. nρCt D4f Tvmax sin δvq

V gz = − Zα αV + g cos γ − ω f cb 2 − δvTc + f a (·) (48)
| {z m }| {zm }
Appl. Sci. 2023, 13, x FOR PEER REVIEW A1 b0 21 of 35
where A1 is the known part of the sink rate control, b0 is the control effectiveness parameter
for the throttle, and fa (·) is the lumped disturbance, which includes the unmodeled dynamics
where ksr is the
and external designed sink
disturbances. rate controller
An NDI control bandwidth.
is designedSimilar to thesink
as a baseline approach in Section
rate controller, as
4.1, to improve the tracking performance, a second-order
shown in Figure 17b. The output of the baseline controller is formulated as follows: ESO is designed to estimate the
disturbance, fa(·), as follows: . .
azE,NDI = −k sr hc − h − A1 (49)
 z 1,sr = β 01,sr ( y − z1,sr ) + z2 ,sr + azE ,d
where ksr is the designed sink rate  control bandwidth. Similar to the approach in Section(50) 4.1,
to improve the tracking performance, ( y − z1,sr )
 z 2 ,sr = β 02a,sr second-order ESO is designed to estimate the
disturbance, fa (·), as follows:
where z1,sr and z2,sr are theestimation
. results of the sink rate and disturbance, and β01,sr and
β02,sr are the ESO’s observation = β 01,sr
z1,srgain, which(y −can be) +
z1,sr z2,sr + azE,das follows [35]:
determined
. (50)
z2,sr = β 02,sr (y − z1,sr )
β01,sr = 2ωo ,sr , β02 ,sr = ωo ,sr 2 (51)
where z1,sr and z2,sr are the estimation results of the sink rate and disturbance, and β01,sr
and β02,sr
where ωo,srare theESO’s
is the ESO’sobservation
observationbandwidth.gain, whichIncan be determined
Equation (50), azE,das
is follows [35]:vertical
the desired
acceleration and can be expressed as follows:
β 01,sr = 2ωo,sr , β 02,sr = ωo,sr 2 (51)
where ω is the ESO’s observation bandwidth. azE ,d = azE ,NDI
In−Equation
z2 ,sr (50), a (52)
is the desired vertical
o,sr zE,d
acceleration
Finally, and can be expressed
the output as follows:
throttle command is as follows:
azE,d = azE,NDI a − z2,sr (52)
δ vTc = zE ,d (53)
Finally, the output throttle command is as follows:
b0
azE,d
δvTc = (53)
b0
4.5. Overall Control Structure
4.5. Overall Control
The overall Structure
control structure of the flying wing UAV, which includes the inner rota-
tionalThe
controller
overall and outer
control translational
structure of thecontroller,
flying wing is UAV,
depicted
whichin Figure 18.the
includes Theinner
methods
rota-
for cross-track
tional controllererror (CTE)translational
and outer control and controller,
vertical offset control in
is depicted areFigure
further
18.discussed in Sec-
The methods for
tion 5.
cross-track error (CTE) control and vertical offset control are further discussed in Section 5.
Figure
Figure 18.
18. Overall
Overall control
control structure for short
structure for short landing.
landing.
5. Short-Landing-Strategy Design
5.1. Short-Landing Procedures
The four-stage short-landing process of the flying wing UAV, shown in Figure 19, is
Appl. Sci. 2023, 13, 3518 21 of 34
5. Short-Landing-Strategy Design
5.1. Short-Landing Procedures
The four-stage short-landing process of the flying wing UAV, shown in Figure 19,
is
Appl. Sci. 2023, 13, x FOR PEER REVIEW designed based on the analyses presented earlier. The process commences at Point
22 of A,
35
where the UAV is at a height of 50 m and a speed of around 64 m/s, with a non-zero sink
rate. The stages of the process are as follows:
Figure
Figure 19. Short-landing process.
19. Short-landing process.
Stage ◦ /s, with the throttle

Stage I:
I: The
The UAVUAV reduces
reduces its
its nozzle
nozzle pitch
pitch angle
angle atat aa rate
rateofof15
15 °/s, with the throttle
fixed
fixed at
at the
the current
current value. The bank
value. The bank angle
angle isis adjusted
adjusted to
to minimize
minimize the the cross-track
cross-track error
error by
by
using the zero sideslip angle command, and the sink rate is controlled by
using the zero sideslip angle command, and the sink rate is controlled by the pitch attitude the pitch attitude
angle, as discussed in Section 4.4. The sink rate command is based on the distance to the
angle, as discussed in Section 4.4. The sink rate command is based on the distance to the
target touchdown point, as explained in Section 5.2. When the nozzle pitch angle reaches
target touchdown point, as explained in Section 5.2. When the nozzle pitch angle reaches
45 degrees, the short landing proceeds to the next stage.
45 degrees, the short landing proceeds to the next stage.
Stage II: In this stage, the lateral and sink rate controls are the same as those in Stage I.
Stage II: In this stage, the lateral and sink rate controls are the same as those in Stage
The nozzle pitch angle is utilized to control the velocity, and the velocity command is set
I. The nozzle pitch angle is utilized to control the velocity, and the velocity command is
to 35 m/s to prevent the inefficiency of the elevator under low dynamic pressure. The
set to 35 m/s to prevent the inefficiency of the elevator under low dynamic pressure. The
short landing advances to Stage III when the velocity decreases to below 36 m/s or the
short landing advances to Stage III when the velocity decreases to below 36 m/s or the
deceleration process lasts more than 20 s.
deceleration process lasts more than 20 s.
Stage III: In this stage, the angle of attack is set to 6 degrees based on the analysis
Stage III: In this stage, the angle of attack is set to 6 degrees based on the analysis
presented in Section 3, and the sink rate is regulated by the throttle, which is different
presented in Section 3, and the sink rate is regulated by the throttle, which is different
from the scenarios in Stages I and II. The lateral control remains the same as in Stage II.
from the scenarios
Additionally, in Stages Ibank
the maximum and angle
II. Thecommand
lateral control remainsto
is restricted the same as below
5 degrees in Stage2 II.
m
Additionally,
above the to
the ground maximum
avoid wingbank angle
strike. Thecommand is restricted
UAV maintains to 5bydegrees
its speed belowangle
nozzle pitch 2m
above the ground
deflection to avoiddown.
until it touches wing strike. The UAV maintains its speed by nozzle pitch angle
deflection until it touches down.
Stage IV: Upon touchdown, the longitudinal control mode switches to pitch angular rate
Stage
control, andIV:
theUpon
commandtouchdown,
is set tothe longitudinal
zero. control
The objective of thismode
stageswitches
is to slowtodown
pitchthe
angular
UAV
rate control, and the command is set to zero. The objective of this stage
by turning off the engine and braking while keeping it on the runway via front-wheel steering.is to slow down
the UAV by turning off the engine and braking while keeping it on the runway via front-
wheel
5.2. steering.
Lateral and Vertical Guidance
The lateral control’s objective is to eliminate the CTE error by using the bank angle
5.2. Lateral
control. Theand Vertical
bank angleGuidance
command is calculated as follows:
The lateral control’s objective is to eliminate
. the CTE
error by using the bank angle
φ = arctan
control. The bank angle commandcis calculateddes ψ V g /g
as follows: (54)

.
(
φc = arctan ψ des Vg g ) (54)
ψdes = ωas ∆ψ
∆ψ = [(ψdes − ψt + π ) mod 2π ] − π (55)
ψ des= =arcsin
ωas Δψ ωy ∆Y/Vg

ψdes
ψ = (ψ des
where ∆Y is the CTE, Vg is theΔground −ψ t + πψ) mod
velocity, 2π  − π angle, and ω is the CTE
(55)
t is the heading y
control bandwidth. The courseψ = arcsin ω ΔY V
angle error,
des
∆ψ,( is
y
limited
g )
between − π and π and is fed
where ΔY is the CTE, Vg is the ground velocity, ψt is the heading angle, and ωy is the
CTE control bandwidth. The course angle error, Δψ, is limited between −π and π and is
fed to a controller with proportional gain, ωas , whose output is the desired rate of change
of course angle.
Appl. Sci. 2023, 13, 3518 22 of 34
To improve the accuracy of the touchdown point, the sink rate command is adjusted
to a controller
online based onwith proportional
the distance fromgain, ω as , whose
the landing point.output is the
The sink ratedesired rate of
command is change
given byof
course angle.
To improve the accuracy of the touchdown point, the sink rate command is adjusted
online based on the distance hc =from
− Vgthe
 ( )
ψ rw −ψ point.
coslanding t 
The( H
 arctan sink ) command is given(56)
Xdesrate by
.
In Figure 20, the projection of the line of sight(LOS) in the direction of the runway is
hc = − Vg cos(|ψrw − ψt |) arctan( H/Xdes ) (56)
denoted as Xdes, the direction angle of the runway is ψrw, and H is the altitude above ground
level (AGL).
In FigureIt is20,
important to noteofthat
the projection thewhen
line ofXdes approaches
sight (LOS) inzero, the sink rate
the direction command
of the runway
tends to infinity
is denoted as Xaccording
des , the to Equation
direction angle (56).
of In
the order
runway to ensure
is ψ rw , the
and landing
H is the gear’s strength,
altitude above
the sink rate
ground levelcommand
(AGL). It is limited
important to atorange
note of −2when
that to 0 m/s,
Xdesand when Xdeszero,
approaches is less
thethan
sink5rate
m,
command
the sink ratetends
commandto infinity according
remains to Equation (56). In order to ensure the landing gear’s
constant.
strength, the sink rate command is limited to a range of −2 to 0 m/s, and when Xdes is less
than 5 m, the sink rate command remains constant.
Figure20.
Figure Sinkrate
20.Sink rateand
andbank
bankangle
anglecommand
commandcalculation.
calculation.
6. Simulation Verification and Discussion

6. Simulation
6.1. Verification
Effectiveness and Discussion
of the Short-Landing Strategy
6.1. Effectiveness
The proposed of the Short-Landing
method’s Strategywas evaluated through numerical simulations,
effectiveness
whichThe proposed
only method’s
considered eventseffectiveness was evaluated
up to the touchdown, through
and the numerical
simulations simulations,
stopped once the
which only considered events up to the touchdown, and the simulations stopped once the
UAV’s landing gear was compressed. Table 4 provides a summary of the control parameters,
UAV’s landing
while Figure 21gear was compressed.
displays Table 4 provides
the three-dimensional a summary
(3D) trajectory of thelanding.
of a short control param-
eters, while Figure 21 displays the three-dimensional (3D) trajectory of a short landing.
Table 4. Inner and outer loop control parameters.
Parameter Explanation Value Parameter Explanation Value

kp p control bandwidth 6 kθ Pitch angle control gain 1.2
kq q control bandwidth 10 kβ Side slip angle control gain 0.8
kr r control bandwidth 6 ksr Sink rate control gain 0.4
ω o,p p observation bandwidth 7 ω o,sr Sink rate observation bandwidth 1
ω o,q q observation bandwidth 5 ωT Velocity control bandwidth 0.25
ω o,r r observation bandwidth 10 ωy CTE control bandwidth 0.15
kφ Roll angle control gain 1.2 ω as Course angle control gain 0.5
The use of thrust vector control in the short-landing process allows for a reduction
in touchdown speed, which, in turn, reduces the required length of the runway. Without
the thrust vector and lift fans, the touchdown speed of the UAV can be calculated [27]
by using the lift coefficient in Figure 6b, resulting in a speed of approximately 46 m/s.
However, with the thrust vector control, the touchdown speed can be reduced to 35
m/s, which is 24% lower than the classic speed. This reduction in touchdown speed
results in a corresponding 24% reduction in the length of the runway, as determined
by the21.
Figure integral relationship
Three-dimensional between
trajectory of velocity and position. As a result, for a given run-
short landing.
Figure 20. Sink rate and bank angle command calculation.
6. Simulation Verification and Discussion

Appl. Sci. 2023, 13, 3518 6.1. Effectiveness of the Short-Landing Strategy 23 of 34
The proposed method’s effectiveness was evaluated through numerical simulations,
which only considered events up to the touchdown, and the simulations stopped once the
UAV’s
way, thelanding gear was
short-landing compressed.
technique Table
allows 4 provides
for the a summary
recovery of the control
of more payload param-
with the same
eters, while Figure 21 displays the three-dimensional
landing speed as the classic landing process. (3D) trajectory of a short landing.
Figure 21. Three-dimensional trajectory of short landing.

Figure 21.
The simulation results are presented in Figure 22, which demonstrates the effectiveness
of the ESO-based angular rate control loops during the short-landing process. Specifically,
the three-axis decoupled control achieved a side slip angle under 1.2 degrees, resulting
in a coordinated turn in Figure 22d,h. The sink rate during the switch to short-landing
mode was approximately −1 m/s, as shown in Figure 22e. During Stage I, the nozzle pitch
angle was increased to 45 degrees, and the throttle command was set to 0.45, as depicted
in Figure 22l,n. As the nozzle angle decreased, the rotation speed increased to 8000 rpm
to balance the moment produced by the vectored thrust. In Stage II, which began when
the nozzle angle reached 45 degrees, the throttle was fixed, and the nozzle pitch angle was
adjusted to reduce the speed. Figure 22l shows that the nozzle pitch angle was increased to
90 degrees in Stage II to reduce the speed from 60 m/s to 35 m/s. During both Stage I and
Stage II, the pitch attitude angle controlled the sink rate, but the sink rate was still affected
by changes in the bank angle, as shown in Figure 22g,e.
Next, the short-landing process entered Stage III, where the nozzle-pitch angle was
used to maintain a constant velocity, the throttle was used to control the sink rate, and a
combination of the elevator and the lift fan controlled the angle of attack. The baseline
angle-of-attack command value was 6 degrees, and Figure 22b shows that there was no
steady-state error in tracking. When the sink-rate control mode was switched to the
ESO-based control in Stage III, the sink-rate controller performed well in terms of its
tracking performance. At around 58 s, a disturbance in the roll angular rate occurred, but it
quickly decayed, as shown in Figure 22e, ultimately improving the landing accuracy. The
maximum deflection of the elevator was approximately −18 degrees, which corresponds to
approximately 60% of its maximum value, as illustrated in Figure 22k.
Appl. Sci. 2023, 13, 3518 24 of 34
(a) (b)
(c) (d)
(e) ( f)
Figure 22. Cont.

Appl. Sci. 2023, 13, 3518 25 of 34
( g) (h)
(i) (j)
(k) (l)
Figure 22. Cont.

Appl.
Appl. Sci.
Sci. 2023,
2023, 13,13, x FOR PEER REVIEW
3518 27ofof3435
26
(m) (n)
(m) (n)
Figure 22. Nominal simulation results of short landing: (a) AGL (Above Ground Level)response, (b)
Figure
Figure Nominal
22.22. Nominalsimulation
simulation results ofshort
results of shortlanding:
landing: (a)(a)
AGLAGL (Above
(Above Ground
Ground Level)response,
Level)response, (b)
angle-of-attack response, (c) velocity response, (d) roll-angle response, (e) sink-rate response, (f)
(b)angle-of-attack
angle-of-attackresponse,
response,(c)(c) velocity
velocity response,
response, (d) (d) roll-angle
roll-angle response,
response, (e) sink-rate
(e) sink-rate response,
response, (f)
roll-angular-rate response, (g) pitch-angle response, (h) side-slip-angle response, (i) pitch-angle-rate
roll-angular-rate
(f)response,
roll-angular-rate response, (g)
response, (g)pitch-angle
pitch-angleresponse, (h) side-slip-angle
response,response, response,
(h) side-slip-angle (i) pitch-angle-rate
response, (i) response,
pitch-angle-
(j) Yaw-angle-rate response, (k) actuator (l) nozzle-pitch-angle (m)
rateresponse,
response, (j) (j)
Yaw-angle-rate
Yaw-angle-rateresponse, (k) actuator
response, (k) response,
actuator (l) nozzle-pitch-angle
response, (l) response,
nozzle-pitch-angle (m)
response,
lift-fan-rotation-speed response, and (n) throttle
lift-fan-rotation-speed response, and (n) throttle response.response.
(m) lift-fan-rotation-speed response, and (n) throttle response.
6.2. Flight Boundary Protectionwith withthe
theAdditional
AdditionalAngle-of-Attack
Angle-of-Attack Command
6.2.6.2. Flight
Flight BoundaryProtection
Boundary Protection with the Additional Angle-of-Attack Command
Command
The
The findings
findingsfrom from the
from the short-landingperformance performance study in in Section 3 indicate thatthat set-
The findings the short-landing
short-landing performance study
study Section
in Section 3 indicate
3 indicate set-
that
ting
ting the
the angle-of-attack
angle-of-attack command
command at
at 66 degrees
degrees cancan achieve
achieve
setting the angle-of-attack command at 6 degrees can achieve good results. However, good
good results.
results. However,
However, when when
certain
whencertain parameters
certainparameters
parameters areare
are perturbed,
perturbed,
perturbed, such
such suchasaslift
aslift coefficient,
coefficient,
lift coefficient,massmass
mass properties,
properties,
properties, and and
and lift
liftlift
fan fan
fan
force,
force, the
the sink
sink rate may
may become
become uncontrollable.
uncontrollable. In In the
the first
first scenario,
scenario,
force, the sink rate may become uncontrollable. In the first scenario, if the lift coefficient if if
the the
lift lift coefficient
coefficient
isisless
isless
less than
thanthanthethe estimated
theestimated
estimated valuevalueor
value or the
orthe UAV’s
theUAV’s
UAV’s weight
weight
weight exceeds
exceeds
exceeds thethethepre-calculated
pre-calculated
pre-calculated value,
value, the the
value,
nominal reference
nominal reference angleangle of
ofattack
attack may
may not
notgenerate
generate enough
enough liftlift
force
forceto maintain
to maintain the the
sink- sink-
the nominal reference angle of attack may not generate enough lift force to maintain the
ratecommand.
rate command. This is is demonstrated
demonstrated ininCase
Case 2 2ininFigure
Figure 23a,
23a,where
where thethe UAV UAV cannot
cannot sus-sus-
sink-rate command. This is demonstrated in Case 2 in Figure 23a, where the UAV cannot
tainits
tain its sink
sink rate and and falls to the ground
groundbefore beforereaching thethetarget point. In the second
sustain its sinkraterate andfallsfalls to
to the
the ground before reaching
reaching target
the target point.
point. InInthethe second
second
scenario, ifif the
scenario, the lift
lift coefficient
coefficient isisgreater
greater than
than the
the CFDCFD calculation
calculation results,
results, thetheUAV UAV moves moves
scenario, if the lift coefficient is greater than the CFD calculation results, the UAV moves
up and down in response to the velocity command, as observed in Case 1. Both of these
upupand anddowndowninin response
response totothethe velocity
velocity command,asas
command, observedinin
observed Case
Case 1. 1.Both
Bothofof these
these
scenarios should
scenarios should be avoided to ensure short-landing precision and safety.
scenarios should bebe avoided
avoided toto ensure
ensure short-landing
short-landing precision
precision andand safety.
safety.
(a) (b)
(a) 23. Simulation and analysis for the angle-of-attack control:
Figure (b)(a) fixed angle-of-attack command
Figure 23. Simulation and analysis for the angle-of-attack control: (a) fixed angle-of-attack command
simulation results, with Case 1 as −20% perturbation of the nominal lift coefficient and Case 2 as
Figure 23.results,
simulation Simulation
with and analysis
Case for
1 as −lift
20% the angle-of-attack
perturbation control: (a) fixed
of the angle-of-attack command
+20% perturbation of the nominal coefficient; and (b) trimnominal liftdifferent
results for coefficient andofCase
angles 2 as
attack
simulation
+20% results, with Case 1 as −20% perturbation of the nominal lift coefficient and Caseand
2 as
andperturbation
sink rates (atof80the nominal
percent liftfan’s
of lift coefficient;
maximum androtation
(b) trimspeed).
results for different angles of attack
+20% perturbation of the nominal lift coefficient; and (b) trim results for different angles of attack
sink rates (at 80 percent of lift fan’s maximum rotation speed).
and sink rates (at 80 percent of lift fan’s maximum rotation speed).
Appl. Sci. 2023, 13, 3518 27 of 34
During short landing, the throttle is used to compensate for the aerodynamic lift force,
while a portion of the lift fan’s control power is used to balance the pitch moment produced
by the vectored thrust. In order to ensure that there is enough pitch controllable margin
to meet the extra pitch control moment requirement, 20% of the lift fan’s control power is
allocated to pitch angular rate control. The trim results, using 80% of the lift fan’s control
power, are presented in Figure 23b.
To prevent the throttle command value from exceeding the lift fan’s trim boundary, an
additional angle-of-attack command is calculated and added to the baseline value. This
additional command can be expressed as follows:
(δvTc − δvTs ) Tvmax sin δvq

∆αc (t) = min max , −∆αc (t − ∆T ) , αcmax − αc (t) (57)
Zα mV
where ∆αc (t) and ∆αc (t − ∆T ) are the current and previous additional angle-of-attack
command values, respectively; αc max is the maximum angle of attack command value,
which is set at 10 degrees; and δvTs is the threshold of the throttle, which is set to 0.6 based
on the trim results in Figure 23b. In Stage III, the current angle-of-attack value is used to
determine the angle-of-attack command value, and the additional angle-of-attack command
value is added to compensate for the perturbations in lift coefficient, mass properties, and
lift fan force, ensuring precise and safe short landings:
αc (t) = α(t) + ∆αc (t) (58)
6.3. Sink-Rate Control under Wind Gust

To evaluate the performance of the ESO-based sink-rate control proposed in Section 4.4,
the implicit NDI controller was chosen for comparison, using the same control bandwidth.
Referring to Equation (41), the control law for sink rate control can be expressed as follows:
Z t . . .
− azE,d = k i hc − h dt − k p h + A1 (59)
0
where kp and ki are the control parameters, which can be determined by referring to
Equations (41) and (42). The wind-gust model is described as follows:

 u g = u2m 1 − cos πx 0 ≤ x ≤ dx
dx (60)
w
w g = m 1 − cos
2
πx
dz 0 ≤ x ≤ dz
where ug is the gust in the u-axis; wg is the gust in the w-axis; x is the integration of velocity;
dx = 120 m and dz = 80 m is the gust length; um = 3.5 m/s; and wm = 3 m/s is the gust
amplitude [40]. As shown in Figure 24, it can be observed that there is no significant
difference in the sink-rate response when there is no wind. However, when the wind is
introduced at 32 s, the ESO-based sink-rate controller exhibits a better disturbance rejection
capability with rapid convergence, while the NDI controller continues to oscillate and
shows a trend of divergence. Therefore, it can be concluded that the ESO-based sink-rate
controller performs better when faced with external disturbances.
Appl. Sci. 2023, 13, 3518 28 of 34
Figure
Figure 24.
24. Sink
Sink rate
ratecontrol
controlunder
underthe
thewind
windgust.
gust.
6.4. Robustness
6.4. Robustness against
againstModel
ModelMismatch
Mismatch
The proposed
proposedshort-landing
short-landingtechnique
techniquemust
mustconsider
considervarious
variousuncertainties
uncertainties and
and mod-
mod-
eling errors
errors to
to ensure
ensureflight
flightsafety.
safety.Examples
Examplesofofthese
theseuncertainties
uncertaintiesinclude
include thrust-induced
thrust-induced
uncertainties[41],
uncertainties [41],ground-effect
ground-effect uncertainties
uncertainties [25],
[25], uncertainties
uncertainties duedue tointeraction
to the the interaction
be-
between lift fans and wings, inertial parameter uncertainties, and uncertainties
tween lift fans and wings, inertial parameter uncertainties, and uncertainties in thrust in thrust
loss
loss induced
induced by the byincrease
the increase in nozzle
in nozzle pitch [42].
pitch angle angleTo [42].
testTo
thetest the robustness
robustness of the pro-
of the proposed
posed technique,
technique, an MC simulation
an MC simulation was conducted
was conducted with 500 simulations,
with 500 simulations, where model where model
parame-
parameters
ters were perturbed
were perturbed in the nonlinear
in the nonlinear model,
model, as shownasinshown
Table in
5. Table 5.
Table 5. Model parameters bias.
Table 5. Model parameters bias.
Parameter
Parameter Explanation
Explanation Bias
BiasRange
Range
M M Mass
Mass ± 1 kg
±1 kg
cgx, cgz
cgx, cgy, cgy, cgz Centre
Centre ofof gravity
gravity ±±5%
5%
qbar Dynamic pressure ±15%
qbar Dynamic pressure ±15%
CLα Lift coefficient derivative ±20%
CLα C Lift coefficient derivative
Roll moment coefficient derivative ±20%
±20%
lβ
CC lβ, C
lδa mδe , C nδr Roll
Roll,moment
pitch, coefficient
and yaw derivative
control derivatives ±±20%
20%
Clδa, Cmδe, CTnδr Roll, pitch, and Engine thrust derivatives
yaw control −20~0%
±20%
T Tf Lift fan
Engine force
thrust −−20~0%
20~0%
Ixx , Iyy , Izz Moment of inertia ±20%
Tf Lift fan force −20~0%
Clp , Cmq Damping moment coefficient derivatives ±50%
Ixx, Iyy, ICzz Moment of inertia
Drag coefficient
±20%
0~20%
D0
Clp, Cmq Damping moment coefficient derivatives ±50%
CD0 the MC simulation, it was
During Drag coefficient
observed 0~20% com-
that, although the angle-of-attack
mand value was set to 6 degrees, the steady-state angle-of-attack response was between
3 and During the MC
9 degrees. simulation,
This shows thatit was
theobserved that, although
flight boundary the angle-of-attack
protection method proposed com-in
mand value was set to 6 degrees, the steady-state angle-of-attack response was
Section 6.2 is effective in protecting the flight boundary. Figure 25a illustrates the steady- between 3
and 9 degrees. This shows that the flight boundary protection method proposed
state angle-of-attack response for the 500 simulations conducted during the MC simulation. in Section
6.2 isTo
effective
preventintheprotecting
throttle the flight boundary.
command value from Figure
going25a illustrates
beyond the boundary,
its trim steady-statean
angle-of-attack response for
additional angle-of-attack the 500 simulations
command is generatedconducted during
by the flight the MC
boundary simulation.
protection method
whenTo prevent
the throttlethe throttle value
command commandexceedsvalue0.6,from goinginbeyond
as shown Figure its trim boundary, an
25b.
additional angle-of-attack
To assess command
the effect of the additional is generated by the
angle-of-attack flight boundary
command protection
on the success rate of
method when the
short landings, throttle
500 command value
MC simulations were exceeds
performed 0.6, without
as shownflight
in Figure 25b. protection
envelope
and using the perturbation range of parameters from Table 5. A landing was considered un-
successful if the distance error between the touchdown point and the intended touchdown
point exceeded 5 m, or if the aircraft failed to land. Among the 500 simulations, 127 short
landings failed due to early touchdown or excessive pitch oscillations during the approach,
as illustrated in Figure 23. The success rate of short landings was found to be 74.6%, and
Figure 26 illustrates the angle-of-attack response of successful landings.
Appl. Sci. 2023, 13, 3518 29 of 34
(a) (b)
Figure 25. States response of the MC simulation: (a) angle-of-attack response and (b) throttle re-
sponse.
To assess the effect of the additional angle-of-attack command on the success rate of
short landings, 500 MC simulations were performed without flight envelope protection
and using the perturbation range of parameters from Table 5. A landing was considered
unsuccessful if the distance error between the touchdown point and the intended touch-
down point exceeded 5 m, or if the aircraft failed to land. Among the 500 simulations, 127
short landings failed due to early touchdown or excessive pitch oscillations during the
(a)
approach, as illustrated in Figure 23. The success rate of short (b)landings was found to be
74.6%, and
Figure Figureresponse
26 illustrates the angle-of-attack response of successful landings.
Figure25.
25.States
States response ofof
thethe
MCMC simulation:
simulation: (a) angle-of-attack
(a) angle-of-attack response response and response.
and (b) throttle (b) throttle re-
sponse.
To assess the effect of the additional angle-of-attack command on the success rate of
short landings, 500 MC simulations were performed without flight envelope protection
and using the perturbation range of parameters from Table 5. A landing was considered
unsuccessful if the distance error between the touchdown point and the intended touch-
down point exceeded 5 m, or if the aircraft failed to land. Among the 500 simulations, 127
short landings failed due to early touchdown or excessive pitch oscillations during the
approach, as illustrated in Figure 23. The success rate of short landings was found to be
74.6%, and Figure 26 illustrates the angle-of-attack response of successful landings.
Figure 26. Angle-of-attack response of successful landing without the flight boundary protection.
Figure 26.
TheAngle-of-attack response
flight trajectories of MC
of the successful landing
simulation without
with the flight
the flight boundary
boundary protection.
protection based
on the additional angle-of-attack command are shown in Figure 27. The UAV follows
The flight
a curved trajectories
trajectory and of the MCasimulation
achieves short landing,with even
the flight
withboundary
parameter protection basedIn
uncertainties.
onaddition,
the additional
there isangle-of-attack
no oscillationcommand
or divergenceare shown
in thein Figure trajectory
landing 27. The UAV follows
during a
the MC
curved trajectory and achieves a short landing, even with parameter uncertainties.
simulation, thus demonstrating the robustness of the proposed control strategy. Moreover, In ad-
dition, there iserrors
the distance no oscillation or were
in all cases divergence
less thanin 4the
m, landing
as shown trajectory
in Figureduring theimprovement
28b. This MC simu-
lation, thus the
enhances demonstrating
likelihood ofthe robustness
successful of the To
landings. proposed
evaluatecontrol strategy.
the landing Moreover,
accuracy, the
the circular
distance errors in all
error probability casesiswere
(CEP) less than
introduced. 4 m,
The as shown
calculated in Figure
result of the 28b.
CEP This improvement
is 0.94 m, as depicted
enhances
in Figure the28a.
likelihood
The mean of successful
and variance landings.
of theTo evaluate
landing the landing
distance accuracy,
error are 0.947 andthe 0.512,
cir-
cular error probability
respectively, as shown (CEP) is introduced.
in Figure The calculated
28b. According to theresult of quality
flying the CEPrequirement
is 0.94 m, as in
depicted
Referencein Figure 28a. The meanlanding
[43], high-precision and variance of the landing distance error are 0.947 and
was achieved.
Figure 26. Angle-of-attack response of successful landing without the flight boundary protection.
The flight trajectories of the MC simulation with the flight boundary protection based
on the additional angle-of-attack command are shown in Figure 27. The UAV follows a
curved trajectory and achieves a short landing, even with parameter uncertainties. In ad-
dition, there is no oscillation or divergence in the landing trajectory during the MC simu-
lation, thus demonstrating the robustness of the proposed control strategy. Moreover, the
distance errors in all cases were less than 4 m, as shown in Figure 28b. This improvement
enhances the likelihood of successful landings. To evaluate the landing accuracy, the cir-
cular error probability (CEP) is introduced. The calculated result of the CEP is 0.94 m, as
depicted in Figure 28a. The mean and variance of the landing distance error are 0.947 and
0.512, respectively, as shown in Figure 28b. According to the flying quality requirement
Appl. Sci. 2023, 13, x FOR PEER REVIEWinReference [43], high-precision landing was achieved. 31 of 35
0.512, respectively, as shown in Figure 28b. According to the flying quality requirement
Appl. Sci. 2023, 13, 3518 30 of 34
in Reference [43], high-precision landing was achieved.
Figure 27. Monte Carlo simulation trajectory.

(a) (b)
Figure28.
Figure 28.Landing
Landingerror
errordiagram:
diagram:(a)
(a)CEP
CEPerror
errorand
and(b)
(b)histogram
histogramwith
witha adistribution
distributionfitfitfor
forCEP.
CEP.
(a) (b)
7.7.Conclusions
Figure Conclusions
28. Landing error diagram: (a) CEP error and (b) histogram with a distribution fit for CEP.
InInsummary,
summary,this
thisstudy
studyinvestigated
investigatedthe
theshort-landing
short-landingtechnology,
technology,and andthe
theproposed
proposed
7. control
Conclusions
control strategy effectively addressed the short-landing problem for a flying-wing UAV.
strategy effectively addressed the short-landing problem for a flying-wing UAV.
The Inresearch
The researchresults
summary, are
areasas
this study
results follows:
investigated
follows: the short-landing technology, and the proposed
control
(1) strategy effectively addressed the short-landing problem for short-landing
a flying-wing UAV.
(1) TheThepaper
paperproposes
proposesa acomplete
completecontrol strategy
control to solve
strategy the
to solve problem
the short-landing for
problem
The research results are as follows:
a flying-wing UAV, using 3BSD thrust-vector control. The frequency-domain-based
for a flying-wing UAV, using 3BSD thrust-vector control. The frequency-domain-
(1) The paperapproach
approach
based proposes ato
is usedis complete
allocate
used control strategy
the liftthe
to allocate fans, to solve
thrust
lift fans, the
vectors,
thrust short-landing
and and
vectors, problemcontrol
aerodynamic
aerodynamic con-
forsurfaces,
a flying-wing UAV,
resulting in ausing 3BSDand
practical thrust-vector control.efficient
computationally The frequency-domain-
solution.
trol surfaces, resulting in a practical and computationally efficient solution.
(2) based
Theapproach is used
equilibrium to allocate
points for thethe lift fans,
velocity and thrust
anglevectors, andare
of attack aerodynamic
determinedcon-
through
troltrim
surfaces, resulting in a practical and computationally efficient solution.
analysis based on the maximum linear and angular acceleration analysis to
establish the command values for velocity and angle of attack.
(3) Different flight modes were considered, and the study designed an ESO-based angular
rate controller and sink rate controller. The three-axis decoupling control is realized
via a SISO design technique, and the ESO-based sink rate controller shows better
disturbance rejection capability when considering gusts of wind.
(4) A flight boundary protection method utilizing an additional angle of attack command
is introduced, which enhances both the precision and safety of landings. Monte Carlo
simulations were conducted to validate the effectiveness of the proposed method.
Appl. Sci. 2023, 13, 3518 31 of 34
8. Future Works
One important research area in short-landing technology is the optimization of the
transition stage when the pitch nozzle begins to increase. This includes exploring methods
such as changing the throttle value when the nozzle pitch angle increases, rather than
keeping a fixed value. In addition, for the potential future application of short-landing
technology to carrier landings, the effects of deck motion and sea state need to be considered
to ensure safe and effective landings.
Author Contributions: Conceptualization, H.L., Y.C. and Z.G.; methodology, H.L. and Y.C.; software,
H.L., Y.Y., Y.B. and Z.G.; validation, H.L., Y.C., Y.B. and Z.G.; formal analysis, H.L., Y.Y., Y.B. and Y.C.;
investigation, H.L. and X.L.; resources, H.L., Y.Y. and Z.G.; data curation, Y.C.; writing—original draft
preparation, H.L. and Y.C.; writing—review and editing, H.L., Y.C., Y.B., Y.Y. and X.L.; visualization, H.L.;
supervision, Y.C. and Z.G. All authors have read and agreed to the published version of the manuscript.
Funding: This work was supported by the National Defense Science and Technology Foundation
Strengthening Program and the Priority Academic Program Development of Jiangsu Higher Educa-
tion Institutions (PAPD).
Institutional Review Board Statement: Ethical review and approval were waived for this study due
to the reason that the study not involving humans or animals.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.
Data Availability Statement: The data used to support the findings of this study are available from
the corresponding author upon reasonable request.
Acknowledgments: The authors would like to deliver their sincere thanks to the editors and anony-
mous reviewers.
Conflicts of Interest: The authors declare that there are no conflicts of interest regarding the publica-
tion of this paper.
Nomenclature
The following nomenclatures are used in this manuscript:
L aerodynamic lift force (N)

D aerodynamic drag force (N)
ϕ roll Euler angle (rad)
θ pitch Euler angle (rad)
V true airspeed (m/s)
p roll angular rate (rad/s)
q pitch angular rate (rad/s)
r yaw angular rate (rad/s)
α angle of attack (rad)
β sideslip angle (rad)
δa the deflection angle of the aileron (rad)
δr the deflection angle of the rudder (rad)
δac commanded deflection angle of the aileron (rad)
δrc commanded deflection angle of the rudder (rad)
δvq 3BSD nozzle pitch angle (rad)
δvr 3BSD nozzle yaw angle (rad)
δvT turbojet engine throttle
vq desired pitch angular acceleration (rad/s2 )
vql low-frequency part of vq (rad/s2 )
vqh high-frequency part of vq (rad/s2 )
vqs the left low-frequency part (rad/s2 )
Appl. Sci. 2023, 13, 3518 32 of 34
ω fcs the command value of the lift fan (rad/s)

ω fcb the baseline rotation speed of the lift fan (rad/s)
ω fcp the parallel allocation part of the lift fan (rad/s)
urde the residual control power of the elevator (rad/s2 )
urfan the residual control power of the lift fan (rad/s2 )
Abbreviations
The following abbreviations are used in this manuscript:
ESO extended state observer

UAV unmanned aerial vehicles
SRVL shipborne rolling vertical landing
ITV integrated thrust-vector
STVL short takeoff and vertical landing
EA eigenstructure assignment
NDI nonlinear dynamic inversion
INDI incremental nonlinear dynamic inversion
ADRC active disturbance rejection control
3BSD three-bearing swivel duct
STOSL short takeoff and short landing
SSD spoiler-slot deflectors
CFD computational fluid dynamics
AAS attainable acceleration set
APAA attainable pitch angular acceleration
BIBO bounded input–bounded output
LFP low-frequency part
HFP high-frequency part
LOS line of sight
AGL above ground level
MC Monte Carlo
CEP circular error probability
SISO single input–single output
References
1. Francis, M.S. Air Vehicle Management with Integrated Thrust-Vector Control. AIAA J. 2018, 56, 4741–4751. [CrossRef]
2. Shanks, G.; Fielding, C.; Andrews, S.; Hyde, R. Flight Control and Handling Research with the VAAC Harrier Aircraft.
Int. J. Control 1994, 59, 291–319. [CrossRef]
3. Postlethwaite, I.; Bates, D.G. Robust Integrated Flight and Propulsion Controller for the Harrier Aircraft. J. Guid. Control. Dyn.
1999, 22, 286–290. [CrossRef]
4. Liang, J.; Bai, J.; Sun, Z.; Liu, S. The Application of Thrust-Vector Control to Damaged Asymmetric Aircraft. In Advances
in Guidance, Navigation and Control; Yan, L., Duan, H., Yu, X., Eds.; Springer: Singapore, 2022; Volume 644, pp. 715–728,
ISBN 9789811581540.
5. Wiegand, C. F-35 Air Vehicle Technology Overview. In Proceedings of the 2018 Aviation Technology, Integration, and Operations
Conference, Atlanta, GA, USA, 25–29 June 2018.
6. Cook, R.; Atkinson, D.; Milla, R.; Revill, N.; Wilson, P. Development of the Shipborne Rolling Vertical Landing (SRVL) Manoeuvre for
the F-35B Aircraft. In Proceedings of the International Powered Lift Conference, Philadelphia, PA, USA, 5–7 October 2010; pp. 5–7.
7. Atkinson, D.C.; Brown, R.; Potts, R.; Bennett, D.; Ward, J.E.; Trott, E. Integration of the F-35 Joint Strike Fighter with the UK
QUEEN ELIZABETH Class Aircraft Carrier. In Proceedings of the 2013 International Powered Lift Conference, Los Angeles, CA,
USA, 12–14 August 2013.
8. F-35 Pilot Makes History with Revolutionary Method of Landing Jets on HMS Queen Elizabeth. Available online: https:
//www.royalnavy.mod.uk/news-and-latest-activity/news (accessed on 8 December 2021).
9. Cheng, Z.; Zhu, J.; Yuan, X. Design of optimal trajectory transition controller for thrust-vectored V/STOL aircraft.
Sci. China Inf. Sci. 2021, 64, 139201. [CrossRef]
10. Gong, Z.; Mao, S.; Wang, Z.; Zhou, Z.; Yang, C.; Li, Z. Control Optimization of Small-Scale Thrust-Vectoring Vertical/Short
Take-Off and Landing Vehicles in Transition Phase. Drones 2022, 6, 129. [CrossRef]
11. Lungu, R.; Lungu, M. Automatic Landing System Using Neural Networks and Radio-Technical Subsystems. Chin. J. Aeronaut.
2017, 30, 399–411. [CrossRef]
Appl. Sci. 2023, 13, 3518 33 of 34
12. Theis, J.; Ossmann, D.; Thielecke, F.; Pfifer, H. Robust Autopilot Design for Landing a Large Civil Aircraft in Crosswind.
Control. Eng. Pract. 2018, 76, 54–64. [CrossRef]
13. Sedlmair, N.; Theis, J.; Thielecke, F. Flight Testing Automatic Landing Control for Unmanned Aircraft Including Curved
Approaches. J. Guid. Control. Dyn. 2022, 45, 726–739. [CrossRef]
14. Ningjun, L.; Zhihao, C.; Yingxun, W.; Jiang, Z. Fast Level-Flight to Hover Mode Transition and Altitude Control in Tiltrotor’s
Landing Operation. Chin. J. Aeronaut. 2021, 34, 181–193.
15. Jiang, J.; Zhong, B.; Fu, S. Influence of Overall Configuration Parameters on Aerodynamic Characteristics of a Blended-Wing-Body
Aircraft. Acta Aeronaut. Astronaut. Sin. 2016, 37, 278.
16. Zhu, Y.; Chen, X.; Li, C. A Moving Frame Trajectory Tracking Method of a Flying-Wing UAV Using the Differential Geometry.
Int. J. Aerosp. Eng. 2016, 2016, 1–9. [CrossRef]
17. Lixin, W.; Zhang, N.; Ting, Y.; Hailiang, L.; Jianghui, Z.; Xiaopeng, J. Three-Axis Coupled Flight Control Law Design for Flying
Wing Aircraft Using Eigenstructure Assignment Method. Chin. J. Aeronaut. 2020, 33, 2510–2526.
18. Walker, G.; Allen, D. X-35B STOVL Flight Control Law Design and Flying Qualities. In Proceedings of the 2002 Biennial
International Powered Lift Conference and Exhibit, Williamsburg, VA, USA, 5–7 November 2002; pp. 1–13.
19. Rui, W.; Zhou, Z.; Yanhang, S. Robust Landing Control and Simulation for Flying Wing UAV. In Proceedings of the 2007 Chinese
Control Conference, Zhangjiajie, China, 26–31 July 2007; pp. 600–604.
20. Zhang, S.; Meng, Q. An Anti-Windup INDI Fault-Tolerant Control Scheme for Flying Wing Aircraft with Actuator Faults.
ISA Trans. 2019, 93, 172–179. [CrossRef]
21. Yin, M.; Chu, Q.P.; Zhang, Y.; Niestroy, M.A.; de Visser, C.C. Probabilistic Flight Envelope Estimation with Application to Unstable
Overactuated Aircraft. J. Guid. Control. Dyn. 2019, 42, 2650–2663. [CrossRef]
22. Lu, K.; Liu, C. Adaptive Control Scheme for UAV Carrier Landing Using Nonlinear Dynamic Inversion. Int. J. Aerosp. Eng. 2019,
2019, 6917393. [CrossRef]
23. Han, J. Active disturbance rejection controller and its applications. Control. Decis. 1998, 13, 19–23.
24. Gao, Z. Active Disturbance Rejection Control: From an Enduring Idea to an Emerging Technology. In Proceedings of the 2015
10th International Workshop on Robot Motion and Control (RoMoCo), Poznań, Poland, 6–8 July 2015; pp. 269–282.
25. Franklin, J.A. Dynamics, Control, and Flying Qualities of V/STOL Aircraft; American Institute of Aeronautics and Astronautics:
Reston, VA, USA, 2002; ISBN 978-1-56347-575-7.
26. Di Francesco, G.; Mattei, M. Modeling and Incremental Nonlinear Dynamic Inversion Control of a Novel Unmanned Tiltrotor.
J. Aircr. 2016, 53, 73–86. [CrossRef]
27. Stevens, B.L.; Lewis, F.L.; Johnson, E.N. Aircraft Control and Simulation: Dynamics, Controls Design, and Autonomous Systems; John
Wiley & Sons: Hoboken, NJ, USA, 2015.
28. Wang, L.; Wang, L.; Jia, Z. Control Features and Application Characteristics of Split Drag Rudder Utilized by Flying Wing. Acta
Aeronaut. Astronaut. Sin. 2011, 32, 1392–1399.
29. Qu, X.; Zhang, W.; Shi, J.; Lyu, Y. A Novel Yaw Control Method for Flying-Wing Aircraft in Low Speed Regime. Aerosp. Sci. Technol.
2017, 69, 636–649. [CrossRef]
30. Franklin, J.A.; Anderson, S.B. V/STOL Maneuverability and Control; National Aeronautics and Space Administration: Washington,
DC, USA, 1984; p. 56.
31. Goman, M.G.; Khramtsovsky, A.V.; Kolesnikov, E.N. Evaluation of Aircraft Performance and Maneuverability by Computation of
Attainable Equilibrium Sets. J. Guid. Control. Dyn. 2008, 31, 329–339. [CrossRef]
32. Durham, W.; Bordignon, K.A.; Beck, R. Aircraft Control Allocation; John Wiley & Sons: Hoboken, NJ, USA, 2017.
33. Li, S.; Yang, J.; Chen, W.-H.; Chen, X. Generalized Extended State Observer Based Control for Systems with Mismatched
Uncertainties. IEEE Trans. Ind. Electron. 2011, 59, 4792–4802. [CrossRef]
34. Zheng, Q.; Gao, Z. Active Disturbance Rejection Control: Between the Formulation in Time and the Understanding in Frequency.
Control. Theory Technol. 2016, 14, 250–259. [CrossRef]
35. Gao, Z. Scaling and Bandwidth-Parameterization Based Controller Tuning. In Proceedings of the 2003 American Control
Conference, Denver, CO, USA, 4–6 June 2003; pp. 4989–4996.
36. Cakiroglu, C.; Van Kampen, E.-J.; Chu, Q.P. Robust Incremental Nonlinear Dynamic Inversion Control Using Angular Ac-
celerometer Feedback. In Proceedings of the 2018 AIAA Guidance, Navigation, and Control Conference, Kissimmee, FL, USA,
8–12 January 2018; p. 1128.
37. Kawaguchi, J.; Miyazawa, Y.; Ninomiya, T. Flight Control Law Design with Hierarchy-Structured Dynamic Inversion Approach. In
Proceedings of the AIAA Guidance, Navigation and Control Conference and Exhibit, Honolulu, HI, USA, 18–21 August 2008; p. 6959.
38. Sieberling, S.; Chu, Q.; Mulder, J. Robust Flight Control Using Incremental Nonlinear Dynamic Inversion and Angular Acceleration
Prediction. J. Guid. Control. Dyn. 2010, 33, 1732–1742. [CrossRef]
39. Lombaerts, T.; Kaneshige, J.; Schuet, S.; Aponso, B.L.; Shish, K.H.; Hardy, G. Dynamic Inversion Based Full Envelope Flight
Control for an EVTOL Vehicle Using a Unified Framework. In Proceedings of the AIAA Scitech 2020 Forum, Orlando, FL, USA,
6–10 January 2020; p. 1619.
40. Lu, Z.; Holzapfel, F. Stability and Performance Analysis for SISO Incremental Flight Control. arXiv 2020, arXiv:2012.00129.
41. Siclari, M.; Migdal, D.; Luzzi Jr, T.; Barche, J.; Palcza, J. Development of Theoretical Models for Jet-Induced Effects on V/STOL
Aircraft. J. Aircr. 1976, 13, 938–944. [CrossRef]
Appl. Sci. 2023, 13, 3518 34 of 34
42. Akturk, A.; Camci, C. Double-Ducted Fan as an Effective Lip Separation Control Concept for Vertical-Takeoff-and-Landing
Vehicles. J. Aircr. 2021, 59, 1–20. [CrossRef]
43. Denison, N.A. Automated Carrier Landing of an Unmanned Combat Aerial Vehicle Using Dynamic Inversion. Master’s Thesis,
Department of the Air Force Air University, Wright-Patterson Air Force Base, Dayton, OH, USA, June 2007; p. 119.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual
author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to
people or property resulting from any ideas, methods, instructions or products referred to in the content.

Applsci 13 03518 v2

Uploaded by

Copyright:

Available Formats

Applsci 13 03518 v2

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Applsci 13 03518 v2

Uploaded by

Copyright:

Available Formats

applied

1 College of Aerospace Engineering, Nanjing University of Aeronautics and Astronautics,

Appl. Sci. 2023, 13, 3518. https://doi.org/10.3390/app13063518 https://www.mdpi.com/journal/applsci

Figure 1. General diagram of the configuration of the subject flyingwing UAV.

Appl. Sci. 2023, 13, x FOR PEER REVIEW 5 of 35

2.3. Thrust Vector Model

Table 3. Effector characteristics of the UAV.

Damping Natural Position

Tvx  −δ vT Tv max  δ vT Tv max cos δ vq cos δ vr 

Figure 4. Longitudinal force analysis.

Figure 5. Pressure contour of the aircraft.

The SSD is used to The SSD is

Appl. Sci. 2023, 13, x FOR PEER REVIEW 11 of 35

3.2. Attainable Linear and Pitch Angular Acceleration

3.2. Attainable Linear and Pitch Angular Acceleration

Equation (13) can be expanded using partial derivatives, as demonstrated in

Figure 14. Complete control diagram for inner and outer

4.2. Control Allocation

otherwise, we have the following:

where we have the following:

 Mδe (δemax − δecs ), vqh < 0

4.3. Angle of Attack Control

where we have the following: 0

α δ vT Tv max sin δ vq nρCtω f D f (45)

Equation (46) can be rewritten as follows:

. nρCt D4f Tvmax sin δvq

Stage ◦ /s, with the throttle

where we have the following:

6. Simulation Verification and Discussion

Parameter Explanation Value Parameter Explanation Value

6. Simulation Verification and Discussion

Figure 21. Three-dimensional trajectory of short landing.

Figure 22. Cont.

Figure 22. Cont.

(δvTc − δvTs ) Tvmax sin δvq

αc (t) = α(t) + ∆αc (t) (58)

6.3. Sink-Rate Control under Wind Gust

Figure 27. Monte Carlo simulation trajectory.

Figure 27. Monte Carlo simulation trajectory.

L aerodynamic lift force (N)

ω fcs the command value of the lift fan (rad/s)

ESO extended state observer

You might also like