# Modeling and Simulation of China C Series Large Aircraft with Microburst

Authors: Jiangyun Wang; Xiao Song; Yongsheng Wang
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

To simulate and analyze a large aircraft with severe wind shear, nonlinear equations for a rigid aircraft under various wind conditions were designed. In this paper, a double-vortex-ring method is adopted to generate a three-dimensional microburst wind field. This model incorporates a flight dynamics module and three flight simulations based on various forms of nonlinear equations. In addition, four primary wind gradients are considered to compute the aerodynamic coefficients. In simulation experiments with a China C series large aircraft, equations using air-relative velocities and inertial velocities are compared and analyzed. Although the nonlinear equations derived using the air-relative velocities should be identical to those using the inertial velocities, the results show that the latter are more accurate in practice because the former bring error into the simulation when estimating the acceleration induced by the wind field.


---

## Body

## 1. Introduction

In the early days of aviation, the influences of a variable wind field upon flight were generally considered to be negligible, even under atmospheric perturbations. Generally, a constant wind velocity was used to estimate the flight performance and plan a flight course. The aviation community did not pay much attention to the threat of variable wind fields until one Boeing-727 flight of Eastern Airlines crashed as a result of serious wind shear in 1975 [1]. Several aviation disasters followed this, all attributed to wind shear, which made the aviation community focus on researching wind shear and its effect on flight. According to a statistical report made by ICAO, from 1968 to 1986, aviation accidents related to weather factors accounted for 30% of the total [2]. Not long before this, the flight accident report published by NTSB showed that, during 1994 to 2003, weather-related flight accidents accounted for 21.3% of the total and half of these were attributed to wind factors [3]. It can be seen that although aeronautical technology is developing rapidly, severe weather, and especially serious wind, has always been a factor in flight accidents. Therefore, the characterization of a wind field with a serious impact on flight has to be further researched, as well as the motion of an aircraft under wind impact.A microburst is one of the most serious wind shear fields because its short duration, small scale, and high intensity cause a series of difficulties in forecasting and detection. It starts with the vertical motion of a column downdraft, which is easily caused by thunderstorm and convective activities. This vertical motion then diverges horizontally close to the ground. The transient variation in wind speed and direction will decrease the performance of an aircraft and result in flight accidents. This fact is acknowledged throughout both the aviation and meteorology communities, especially during takeoff and approach stages [4]. Since the late 1970s, a series of projects such as NIMROD and JAWS have been set up to research the characterization and hazards of wind shear using actually measured data and digital flight data records [1, 5–7]. As an extension and intension of these studies, mathematical models for simulating microburst have been modified to determine the influence of each parameter extracted from wind characterization data on the flight performance. Typically, models based on fluid dynamics have been developed because of their simplicity and high precision. Among these, the doublet sheet model [8], vortex-ring model, and vortex section superposition model [9, 10] are usually adopted to construct real-time wind fields. In fact, the first two are identical in nature and have the same results, although the doublet of the first model is perpendicular to the vortex pair in the vortex-ring model. Moreover, the doublet sheet model is usually fit to theoretical research because of the heavy burden of computing its integral operations. However, the third one is merely an approximation of the vortex-ring model. In summary, a double-vortex-ring algorithm based on the vortex-ring model, according to [4], will be chosen to establish the wind field in this paper.Regarding the motion equations for an aircraft under a wind field, Etkin [11] first considered a constant wind velocity and derived equations in different coordinates. However, the deviation of an aircraft’s motion under a constant wind field can easily be modified, which is impossible under variable wind. Thus, Frost and Bowles [12] derived equations in detail under a variable wind field using air-relative velocities. Compared to those without disturbance, their equations incorporated the terms of the wind velocity vector and its derivative. They also compared the different forms of wind terms between different coordinates. In addition, they applied their equations to compute the response of an aircraft under wind and establish flight simulators with a wind effect module [1]. It is remarkable that the absolute derivative of the wind velocity vector, which was called the effective force multiplied by the mass of an aircraft in their thesis, was attributed to using air-relative velocities to derive the force equations. It is identical to the acceleration induced by the relative motion of the air flow to the inertial coordinates.Moreover, a wind field can have a serious impact on flight dynamics when the acceleration induced by wind movement has the same order of magnitude as the gravitational acceleration [12]. Yuan et al. [13] presented a tracker to minimize the effect of wind on the impact-point prediction (IPP) performance. There is no doubt that all of this work laid the foundation for research on the motion characterization of large vehicles under a wind field. Furthermore, Rauw derived the force equations, under a wind field, for the airspeed, incidence, and sideslip angle according to the relationship between the airspeed and air-relative velocity vector [14]. Some scholars recommend this kind of force equation, especially when variables like the derivative of the incidence are included in the aerodynamic model.In conclusion, many researchers prefer the equations derived using air-relative velocities for two reasons: the aerodynamic forces are defined relative to the air flow and it is convenient to explain the effect of wind velocity intuitively. Actually, the equations derived using the inertial velocities incorporate the impact of the wind velocity when computing the aerodynamic coefficients. The changes in the aerodynamic forces are equal to the forces generated by the wind velocity, and the changes in the aerodynamic moments are equivalent to the rotational motion induced by wind gradients. In theory, both forms of equations are the same since they both obey Newton’s second law. However, the practical consequences are usually unpredictable. In [1], the author regarded the kind of equations to be chosen as an individual preference.Therefore, the different forms of the motion equations for a rigid aircraft under a wind field are first reviewed in this paper. Next, we explain the influence of wind gradients upon the angular rates and motion of an aircraft. Then, a microburst wind field is built according to [4, 15] and included in setting up a simulation platform. The results of experiments are used to explain the differences between the wind influence mechanisms in the two forms of equations.


## 2. Modeling of Large Aircraft under Wind Field

### 2.1. Motion Equations of Large Aircraft under Wind Field

In this paper, the North-East-Down axis system is regarded as an inertial reference frame, and all the equations are derived using body coordinates. In addition, six degree of freedom motion equations under a wind field are developed on the basis of the following four hypotheses:First, the earth is flat and has no rotation.
Second, the aircraft is regarded as a rigid body and its mass is constant.
Third, the wind field is called a “Taylor Frozen Field,” so it is constant with respect to time.
Last, the aircraft body, including both its geometric profile and built-in mass distribution, is symmetrical about theX
O
Z plane, so I
x
y
=
I
y
z
=
0.


#### 2.1.1. Nonlinear Equations without Perturbations

Typically, the forces acting on an aircraft primarily include the liftL, drag D, side force C, gravity G, and thrust T, which have the directions shown in Figure 1. Among these, the lift, drag, and side forces are called aerodynamic forces because they are caused by the motion of the aircraft relative to the air flow. Thus, their directions are defined in accordance with the direction of the air-relative velocity vector and generally computed in the air coordinates frame. In general, thrust is regarded as acting in the symmetric plane of the rigid aircraft with an inclination denoted by σ.Figure 1
Diagram of forces and velocities.LetV
→
=
u
v
w
T denote the air-relative velocity vector at the center of gravity, and let Ω
→
=
p
q
r
T be the angular velocity vector. Without perturbations, the inertial velocity is identical to the air-relative velocity. Then, the total force acting on the aircraft is represented in body coordinates as follows:(1)
F
→
=
F
x
F
y
F
z
=
T
cos
⁡
σ
0
T
sin
⁡
σ
+
C
B
A
-
D
C
-
L
+
C
B
E
0
0
m
g
,where the symbols C
B
A and C
B
E are the transform matrices from air coordinates and inertial coordinates to body coordinates, respectively. The minus sign in the drag and lift components means that their directions are the opposite to those of the coordinates adopted. According to Newton’s second law, the dynamic equation of force in body coordinates is as follows:(2)
m
V
→
˙
+
Ω
→
×
V
→
=
F
→
.Its components are rewritten as follows:(3)
u
˙
=
F
x
m
+
r
v
-
q
w
v
˙
=
F
y
m
+
p
w
-
r
u
w
˙
=
F
z
m
+
q
u
-
p
v
.LetV
as, α, β denote the airspeed, incidence, and sideslip angle, respectively, which correspond to the air-relative velocity by the following relationships:(4)
V
as
=
u
2
+
v
2
+
w
2
,
α
=
arctan
⁡
w
u
,
β
=
arcsin
⁡
v
V
as
.By differentiating the above with respect to time and replacingu
˙, v
˙, w
˙ with (3), we can obtain the equations of force in [14] as follows:(5)
V
˙
as
=
1
m
F
x
cos
⁡
α
cos
⁡
β
+
F
y
sin
⁡
β
+
F
z
sin
⁡
α
cos
⁡
β
,
α
˙
=
-
F
x
sin
⁡
α
+
F
z
cos
⁡
α
m
V
as
cos
⁡
β
+
q
-
tan
⁡
β
p
cos
⁡
α
+
r
sin
⁡
α
,
β
˙
=
1
m
V
as
-
F
x
cos
⁡
α
sin
⁡
β
+
F
y
cos
⁡
β
-
F
z
sin
⁡
α
sin
⁡
β
+
p
sin
⁡
α
-
r
cos
⁡
α
.Here, we can state that (5) is exactly the same as (3), in both theory and practice, as long as the forces are equal.Second, the equations of navigation are as follows:(6)
x
˙
=
u
cos
⁡
θ
cos
⁡
ψ
+
v
sin
⁡
ϕ
sin
⁡
θ
cos
⁡
ψ
-
cos
⁡
ϕ
sin
⁡
ψ
+
w
cos
⁡
ϕ
sin
⁡
θ
cos
⁡
ψ
+
sin
⁡
ϕ
sin
⁡
ψ
,
y
˙
=
u
cos
⁡
θ
sin
⁡
ψ
+
v
sin
⁡
ϕ
sin
⁡
θ
sin
⁡
ψ
+
cos
⁡
ϕ
cos
⁡
ψ
+
w
cos
⁡
ϕ
sin
⁡
θ
sin
⁡
ψ
-
sin
⁡
ϕ
cos
⁡
ψ
,
z
˙
=
-
u
sin
⁡
θ
+
v
sin
⁡
ψ
cos
⁡
θ
+
w
cos
⁡
ψ
cos
⁡
θ
.Third, the dynamic equations of moment in body coordinates are directly given without deduction:(7)
p
˙
=
1
I
x
x
I
z
z
-
I
x
z
2
I
x
x
-
I
y
y
+
I
z
z
I
x
z
p
q
+
I
y
y
I
z
z
-
I
z
z
2
-
I
x
z
2
q
r
+
I
z
z
l
+
I
x
z
N
,
q
˙
=
1
I
y
y
I
z
z
-
I
x
x
p
r
+
I
x
z
r
2
-
p
2
+
M
,
r
˙
=
1
I
x
x
I
z
z
-
I
x
z
2
I
x
x
2
-
I
x
x
I
y
y
+
I
x
z
2
p
q
+
-
I
x
x
+
I
y
y
-
I
z
z
I
x
z
q
r
+
I
x
z
l
+
I
x
x
N
.In the above,l, M, N are the components of the total external torque vector in body coordinates. I
x
x, I
y
y, I
z
z denote the moments of inertia, and I
x
z is the product of inertia of the rigid aircraft.Finally, the kinematic equations of rotation are as follows:(8)
ψ
˙
=
q
sin
⁡
ϕ
+
r
cos
⁡
ϕ
cos
⁡
θ
,
θ
˙
=
q
cos
⁡
ϕ
-
r
sin
⁡
ϕ
,
ϕ
˙
=
p
+
tan
⁡
θ
·
q
sin
⁡
ϕ
+
r
cos
⁡
ϕ
.


#### 2.1.2. Equations Derived Using Air-Relative Velocities under Wind Effect

A vector triangle, namely,V
→
e
=
V
→
+
W
→, is the most fundamental relationship between velocities under wind. In order to distinguish these vectors, we use the subscript e to denote the vectors relative to the inertial coordinates and subscript w to denote the variables related to the wind velocity. Thus, let u
e, v
e, w
e denote the components of inertial velocity vector V
→
e in body coordinates, and let W
→
=
u
w
v
w
w
w
T denote the wind vector. To set up the nonlinear equations using air-relative velocities, we can substitute V
→
e in the relationship V
→
e
=
V
→
+
W
→ for the term V
→ in (2), as follows:(9)
m
V
→
˙
+
W
→
˙
+
Ω
→
×
V
→
+
W
→
=
F
→
.After arranging, we have the following:(10)
V
→
˙
=
F
→
m
-
Ω
→
×
V
→
-
W
→
˙
+
Ω
→
×
W
→
,where the term W
→
˙
+
Ω
→
×
W
→ is the absolute derivative of the wind velocity vector and called the effective force to be multiplied by the mass in [1, 12].The components of the above are as follows:(11)
u
˙
=
F
x
m
+
r
v
-
q
w
-
u
˙
w
+
r
v
w
-
q
w
w
,
v
˙
=
F
y
m
+
p
w
-
r
u
-
v
˙
w
+
p
w
w
-
r
u
w
,
w
˙
=
F
z
m
+
q
u
-
p
v
-
w
˙
w
+
q
u
w
-
p
v
w
.We can directly obtain the airspeed and incidence from the above through the mutual relationship in (4). The derivatives of the wind are as follows:(12)
u
˙
w
=
∂
u
w
∂
t
+
∂
u
w
∂
x
x
˙
e
+
∂
u
w
∂
y
y
˙
e
+
∂
u
w
∂
z
z
˙
e
,
v
˙
w
=
∂
v
w
∂
t
+
∂
v
w
∂
x
x
˙
e
+
∂
v
w
∂
y
y
˙
e
+
∂
v
w
∂
z
z
˙
e
,
w
˙
w
=
∂
w
w
∂
t
+
∂
w
w
∂
x
x
˙
e
+
∂
w
w
∂
y
y
˙
e
+
∂
w
w
∂
z
z
˙
e
.Since the motion of the air flow is slow enough relative to the airspeed of an aircraft, the wind velocity is generally regarded as steady with respect to time; namely,∂
u
w
/
∂
t
=
∂
v
w
/
∂
t
=
∂
w
w
/
∂
t
=
0. This is why we assume that the wind field satisfies the hypothesis of Taylor’s Frozen Field. The residual terms all include the wind gradients, as will be addressed later.Let(13)
F
→
w
=
F
w
x
F
w
y
F
w
z
=
m
u
˙
w
+
r
v
w
-
q
w
w
m
v
˙
w
+
p
w
w
-
r
u
w
m
w
˙
w
+
q
u
w
-
p
v
wdenote the effective force generated by the wind field. Then, the state equations for the airspeed, incidence, and sideslip angle in [14] are merely cited as follows:(14)
V
˙
as
=
1
m
F
x
-
F
w
x
cos
⁡
α
cos
⁡
β
+
F
y
-
F
w
y
sin
⁡
β
+
F
z
-
F
w
z
sin
⁡
α
cos
⁡
β
,
α
˙
=
-
F
x
-
F
w
x
sin
⁡
α
+
F
z
-
F
w
z
cos
⁡
α
m
V
as
cos
⁡
β
+
q
-
tan
⁡
β
p
cos
⁡
α
+
r
sin
⁡
α
,
β
˙
=
1
m
V
as
-
F
x
-
F
w
x
cos
⁡
α
sin
⁡
β
+
F
y
-
F
w
y
cos
⁡
β
-
F
z
-
F
w
z
sin
⁡
α
sin
⁡
β
+
p
sin
⁡
α
-
r
cos
⁡
α
.Compared with (5), the dynamic equations of force above are generally the same in form, only to subtract the effective force generated by the wind acceleration.Equations of navigation are obtained usingV
→
+
W
→ to replace V
→ in (6), as follows:(15)
x
˙
e
=
u
+
u
w
cos
⁡
θ
cos
⁡
ψ
+
v
+
v
w
sin
⁡
ϕ
sin
⁡
θ
cos
⁡
ψ
-
cos
⁡
ϕ
sin
⁡
ψ
+
w
+
w
w
cos
⁡
ϕ
sin
⁡
θ
cos
⁡
ψ
+
sin
⁡
ϕ
sin
⁡
ψ
,
y
˙
e
=
u
+
u
w
cos
⁡
θ
sin
⁡
ψ
+
v
+
v
w
sin
⁡
ϕ
sin
⁡
θ
sin
⁡
ψ
+
cos
⁡
ϕ
cos
⁡
ψ
+
w
+
w
w
cos
⁡
ϕ
sin
⁡
θ
sin
⁡
ψ
-
sin
⁡
ϕ
cos
⁡
ψ
,
z
˙
e
=
-
u
+
u
w
sin
⁡
θ
+
v
+
v
w
sin
⁡
ψ
cos
⁡
θ
+
w
+
w
w
cos
⁡
ψ
cos
⁡
θ
.Since the velocities are not explicit in the dynamic and kinematic equations of rotation, the effect of the wind field on these equations is transitive through the aerodynamic moment coefficients. Even under wind perturbations, the equations are all the same as those without disturbance. Thus, it is not necessary to rewrite them again here.


#### 2.1.3. Equations Derived Using Inertial Velocities under Wind Effect

All of the nonlinear equations under a wind condition derived using inertial velocities have the same form as those without perturbations, withV
→
e simply substituted for V
→. Thus, we can easily obtain the dynamic equations of force using the inertial velocities:(16)
u
˙
e
=
F
x
m
+
r
v
e
-
q
w
e
,
v
˙
e
=
F
y
m
+
p
w
e
-
r
u
e
,
w
˙
e
=
F
z
m
+
q
u
e
-
p
v
e
.Then, the equations of navigation are as follows:(17)
x
˙
e
=
u
e
cos
⁡
θ
cos
⁡
ψ
+
v
e
sin
⁡
ϕ
sin
⁡
θ
cos
⁡
ψ
-
cos
⁡
ϕ
sin
⁡
ψ
+
w
e
cos
⁡
ϕ
sin
⁡
θ
cos
⁡
ψ
+
sin
⁡
ϕ
sin
⁡
ψ
,
y
˙
e
=
u
e
cos
⁡
θ
sin
⁡
ψ
+
v
e
sin
⁡
ϕ
sin
⁡
θ
sin
⁡
ψ
+
cos
⁡
ϕ
cos
⁡
ψ
+
w
e
cos
⁡
ϕ
sin
⁡
θ
sin
⁡
ψ
-
sin
⁡
ϕ
cos
⁡
ψ
,
z
˙
e
=
-
u
e
sin
⁡
θ
+
v
e
sin
⁡
ψ
cos
⁡
θ
+
w
e
cos
⁡
ψ
cos
⁡
θ
.In accordance with the velocity triangle relationship, the components of the air-relative velocity vector in body coordinates are as follows:(18)
u
v
w
=
u
e
v
e
w
e
-
u
w
v
w
w
w
.Consequently, the airspeed and incidence can be obtained as follows:(19)
V
as
=
u
2
+
v
2
+
w
2
,
α
=
arctan
⁡
w
u
,
β
=
arcsin
⁡
v
V
as
.Differentiating (19), we can obtain the state equations of the flight parameters and find them the same as (14).


### 2.2. Modeling of Aerodynamic Coefficients

Now, the data for modeling the aerodynamic coefficients can primarily be obtained from flight test records, wind tunnel experiment data, and aided computation methods. The nondimensional aerodynamic coefficients are typically dependent on such variables such as the incidence, sideslip angle, elevator deflection, rotational angular rates, rudder deflection, and aileron deflection. Of these, each aerodynamic derivative is given by a one- or two-dimensional table with respect to the current height and Mach and computed through linear interpolation when given the flight state parameters as inputs. The coefficients model adopted here uses components in body coordinates:(20)
C
x
=
C
x
0
+
C
x
α
α
+
C
x
δ
r
δ
r
,
C
y
=
C
y
β
β
+
C
y
α
β
α
β
+
C
y
δ
r
δ
r
+
C
y
δ
a
δ
a
,
C
z
=
C
z
α
α
+
C
z
δ
e
δ
e
,
C
l
=
C
l
β
β
+
C
l
β
α
α
β
+
C
l
p
p
b
2
V
+
C
l
r
r
b
2
V
+
C
l
δ
r
δ
r
+
C
l
δ
a
δ
a
+
C
l
δ
a
α
δ
a
α
,
C
m
=
C
m
0
+
C
m
α
α
+
C
m
q
q
c
-
V
+
C
m
δ
e
δ
e
+
C
m
r
r
b
2
V
,
C
n
=
C
n
0
+
C
n
β
α
α
β
+
C
n
β
β
+
C
n
p
p
b
2
V
+
C
n
r
r
b
2
V
+
C
n
δ
r
δ
r
+
C
n
δ
a
δ
a
.


### 2.3. Modeling of Wind Gradients

Since different parts of a rigid aircraft under a variable wind field may be influenced by different wind velocity vectors simultaneously, this dissimilarity would result in changes in the aerodynamic forces and moments acting on the aircraft [1]. Evans et al. [16] summarized the influence of spatial wind gradients upon the aerodynamic coefficients. In their thesis, a modified vortex-lattice method that took wind gradients into consideration was adopted to compute the aerodynamic coefficients of a wing and horizontal tail. Through their analyses and comparisons, they concluded that, under a microburst wind shear, wind gradients have such a significant influence upon the aerodynamic moments that they cannot be negligible. However, the impact on the aerodynamic forces is relatively unimportant. Actually, there are nine terms for wind gradients in spatial wind fields, and it is difficult to determine which is solely responsible for pitching, rolling, or yawing because of the interaction between these motions. In addition, the majority of them have no significant impact on the motion of an aircraft, especially under a small angle of attack. Etkin pointed out that a nonuniform distribution of air flow along the fuselage and wingspan would cause additional aerodynamic forces and moments. For the sake of simplification, based on the assumption that the wind velocity distribution is uniform along the coordinates, he gave four primary wind gradients, namely, ∂
w
w
/
∂
y, ∂
w
w
/
∂
x, ∂
v
w
/
∂
x, and ∂
u
w
/
∂
y. These will cause negative rolling, positive pitching, and negative and positive yawing, respectively, in body coordinates [17]. Thus, the rotation angular velocity used to compute the aerodynamic coefficients can be modified as follows:(21)
p
eff
q
eff
r
eff
=
p
q
r
+
-
∂
w
w
∂
y
∂
w
w
∂
x
∂
u
w
∂
y
-
∂
v
w
∂
x
.


## 2.1. Motion Equations of Large Aircraft under Wind Field

In this paper, the North-East-Down axis system is regarded as an inertial reference frame, and all the equations are derived using body coordinates. In addition, six degree of freedom motion equations under a wind field are developed on the basis of the following four hypotheses:First, the earth is flat and has no rotation.
Second, the aircraft is regarded as a rigid body and its mass is constant.
Third, the wind field is called a “Taylor Frozen Field,” so it is constant with respect to time.
Last, the aircraft body, including both its geometric profile and built-in mass distribution, is symmetrical about theX
O
Z plane, so I
x
y
=
I
y
z
=
0.


### 2.1.1. Nonlinear Equations without Perturbations

Typically, the forces acting on an aircraft primarily include the liftL, drag D, side force C, gravity G, and thrust T, which have the directions shown in Figure 1. Among these, the lift, drag, and side forces are called aerodynamic forces because they are caused by the motion of the aircraft relative to the air flow. Thus, their directions are defined in accordance with the direction of the air-relative velocity vector and generally computed in the air coordinates frame. In general, thrust is regarded as acting in the symmetric plane of the rigid aircraft with an inclination denoted by σ.Figure 1
Diagram of forces and velocities.LetV
→
=
u
v
w
T denote the air-relative velocity vector at the center of gravity, and let Ω
→
=
p
q
r
T be the angular velocity vector. Without perturbations, the inertial velocity is identical to the air-relative velocity. Then, the total force acting on the aircraft is represented in body coordinates as follows:(1)
F
→
=
F
x
F
y
F
z
=
T
cos
⁡
σ
0
T
sin
⁡
σ
+
C
B
A
-
D
C
-
L
+
C
B
E
0
0
m
g
,where the symbols C
B
A and C
B
E are the transform matrices from air coordinates and inertial coordinates to body coordinates, respectively. The minus sign in the drag and lift components means that their directions are the opposite to those of the coordinates adopted. According to Newton’s second law, the dynamic equation of force in body coordinates is as follows:(2)
m
V
→
˙
+
Ω
→
×
V
→
=
F
→
.Its components are rewritten as follows:(3)
u
˙
=
F
x
m
+
r
v
-
q
w
v
˙
=
F
y
m
+
p
w
-
r
u
w
˙
=
F
z
m
+
q
u
-
p
v
.LetV
as, α, β denote the airspeed, incidence, and sideslip angle, respectively, which correspond to the air-relative velocity by the following relationships:(4)
V
as
=
u
2
+
v
2
+
w
2
,
α
=
arctan
⁡
w
u
,
β
=
arcsin
⁡
v
V
as
.By differentiating the above with respect to time and replacingu
˙, v
˙, w
˙ with (3), we can obtain the equations of force in [14] as follows:(5)
V
˙
as
=
1
m
F
x
cos
⁡
α
cos
⁡
β
+
F
y
sin
⁡
β
+
F
z
sin
⁡
α
cos
⁡
β
,
α
˙
=
-
F
x
sin
⁡
α
+
F
z
cos
⁡
α
m
V
as
cos
⁡
β
+
q
-
tan
⁡
β
p
cos
⁡
α
+
r
sin
⁡
α
,
β
˙
=
1
m
V
as
-
F
x
cos
⁡
α
sin
⁡
β
+
F
y
cos
⁡
β
-
F
z
sin
⁡
α
sin
⁡
β
+
p
sin
⁡
α
-
r
cos
⁡
α
.Here, we can state that (5) is exactly the same as (3), in both theory and practice, as long as the forces are equal.Second, the equations of navigation are as follows:(6)
x
˙
=
u
cos
⁡
θ
cos
⁡
ψ
+
v
sin
⁡
ϕ
sin
⁡
θ
cos
⁡
ψ
-
cos
⁡
ϕ
sin
⁡
ψ
+
w
cos
⁡
ϕ
sin
⁡
θ
cos
⁡
ψ
+
sin
⁡
ϕ
sin
⁡
ψ
,
y
˙
=
u
cos
⁡
θ
sin
⁡
ψ
+
v
sin
⁡
ϕ
sin
⁡
θ
sin
⁡
ψ
+
cos
⁡
ϕ
cos
⁡
ψ
+
w
cos
⁡
ϕ
sin
⁡
θ
sin
⁡
ψ
-
sin
⁡
ϕ
cos
⁡
ψ
,
z
˙
=
-
u
sin
⁡
θ
+
v
sin
⁡
ψ
cos
⁡
θ
+
w
cos
⁡
ψ
cos
⁡
θ
.Third, the dynamic equations of moment in body coordinates are directly given without deduction:(7)
p
˙
=
1
I
x
x
I
z
z
-
I
x
z
2
I
x
x
-
I
y
y
+
I
z
z
I
x
z
p
q
+
I
y
y
I
z
z
-
I
z
z
2
-
I
x
z
2
q
r
+
I
z
z
l
+
I
x
z
N
,
q
˙
=
1
I
y
y
I
z
z
-
I
x
x
p
r
+
I
x
z
r
2
-
p
2
+
M
,
r
˙
=
1
I
x
x
I
z
z
-
I
x
z
2
I
x
x
2
-
I
x
x
I
y
y
+
I
x
z
2
p
q
+
-
I
x
x
+
I
y
y
-
I
z
z
I
x
z
q
r
+
I
x
z
l
+
I
x
x
N
.In the above,l, M, N are the components of the total external torque vector in body coordinates. I
x
x, I
y
y, I
z
z denote the moments of inertia, and I
x
z is the product of inertia of the rigid aircraft.Finally, the kinematic equations of rotation are as follows:(8)
ψ
˙
=
q
sin
⁡
ϕ
+
r
cos
⁡
ϕ
cos
⁡
θ
,
θ
˙
=
q
cos
⁡
ϕ
-
r
sin
⁡
ϕ
,
ϕ
˙
=
p
+
tan
⁡
θ
·
q
sin
⁡
ϕ
+
r
cos
⁡
ϕ
.


### 2.1.2. Equations Derived Using Air-Relative Velocities under Wind Effect

A vector triangle, namely,V
→
e
=
V
→
+
W
→, is the most fundamental relationship between velocities under wind. In order to distinguish these vectors, we use the subscript e to denote the vectors relative to the inertial coordinates and subscript w to denote the variables related to the wind velocity. Thus, let u
e, v
e, w
e denote the components of inertial velocity vector V
→
e in body coordinates, and let W
→
=
u
w
v
w
w
w
T denote the wind vector. To set up the nonlinear equations using air-relative velocities, we can substitute V
→
e in the relationship V
→
e
=
V
→
+
W
→ for the term V
→ in (2), as follows:(9)
m
V
→
˙
+
W
→
˙
+
Ω
→
×
V
→
+
W
→
=
F
→
.After arranging, we have the following:(10)
V
→
˙
=
F
→
m
-
Ω
→
×
V
→
-
W
→
˙
+
Ω
→
×
W
→
,where the term W
→
˙
+
Ω
→
×
W
→ is the absolute derivative of the wind velocity vector and called the effective force to be multiplied by the mass in [1, 12].The components of the above are as follows:(11)
u
˙
=
F
x
m
+
r
v
-
q
w
-
u
˙
w
+
r
v
w
-
q
w
w
,
v
˙
=
F
y
m
+
p
w
-
r
u
-
v
˙
w
+
p
w
w
-
r
u
w
,
w
˙
=
F
z
m
+
q
u
-
p
v
-
w
˙
w
+
q
u
w
-
p
v
w
.We can directly obtain the airspeed and incidence from the above through the mutual relationship in (4). The derivatives of the wind are as follows:(12)
u
˙
w
=
∂
u
w
∂
t
+
∂
u
w
∂
x
x
˙
e
+
∂
u
w
∂
y
y
˙
e
+
∂
u
w
∂
z
z
˙
e
,
v
˙
w
=
∂
v
w
∂
t
+
∂
v
w
∂
x
x
˙
e
+
∂
v
w
∂
y
y
˙
e
+
∂
v
w
∂
z
z
˙
e
,
w
˙
w
=
∂
w
w
∂
t
+
∂
w
w
∂
x
x
˙
e
+
∂
w
w
∂
y
y
˙
e
+
∂
w
w
∂
z
z
˙
e
.Since the motion of the air flow is slow enough relative to the airspeed of an aircraft, the wind velocity is generally regarded as steady with respect to time; namely,∂
u
w
/
∂
t
=
∂
v
w
/
∂
t
=
∂
w
w
/
∂
t
=
0. This is why we assume that the wind field satisfies the hypothesis of Taylor’s Frozen Field. The residual terms all include the wind gradients, as will be addressed later.Let(13)
F
→
w
=
F
w
x
F
w
y
F
w
z
=
m
u
˙
w
+
r
v
w
-
q
w
w
m
v
˙
w
+
p
w
w
-
r
u
w
m
w
˙
w
+
q
u
w
-
p
v
wdenote the effective force generated by the wind field. Then, the state equations for the airspeed, incidence, and sideslip angle in [14] are merely cited as follows:(14)
V
˙
as
=
1
m
F
x
-
F
w
x
cos
⁡
α
cos
⁡
β
+
F
y
-
F
w
y
sin
⁡
β
+
F
z
-
F
w
z
sin
⁡
α
cos
⁡
β
,
α
˙
=
-
F
x
-
F
w
x
sin
⁡
α
+
F
z
-
F
w
z
cos
⁡
α
m
V
as
cos
⁡
β
+
q
-
tan
⁡
β
p
cos
⁡
α
+
r
sin
⁡
α
,
β
˙
=
1
m
V
as
-
F
x
-
F
w
x
cos
⁡
α
sin
⁡
β
+
F
y
-
F
w
y
cos
⁡
β
-
F
z
-
F
w
z
sin
⁡
α
sin
⁡
β
+
p
sin
⁡
α
-
r
cos
⁡
α
.Compared with (5), the dynamic equations of force above are generally the same in form, only to subtract the effective force generated by the wind acceleration.Equations of navigation are obtained usingV
→
+
W
→ to replace V
→ in (6), as follows:(15)
x
˙
e
=
u
+
u
w
cos
⁡
θ
cos
⁡
ψ
+
v
+
v
w
sin
⁡
ϕ
sin
⁡
θ
cos
⁡
ψ
-
cos
⁡
ϕ
sin
⁡
ψ
+
w
+
w
w
cos
⁡
ϕ
sin
⁡
θ
cos
⁡
ψ
+
sin
⁡
ϕ
sin
⁡
ψ
,
y
˙
e
=
u
+
u
w
cos
⁡
θ
sin
⁡
ψ
+
v
+
v
w
sin
⁡
ϕ
sin
⁡
θ
sin
⁡
ψ
+
cos
⁡
ϕ
cos
⁡
ψ
+
w
+
w
w
cos
⁡
ϕ
sin
⁡
θ
sin
⁡
ψ
-
sin
⁡
ϕ
cos
⁡
ψ
,
z
˙
e
=
-
u
+
u
w
sin
⁡
θ
+
v
+
v
w
sin
⁡
ψ
cos
⁡
θ
+
w
+
w
w
cos
⁡
ψ
cos
⁡
θ
.Since the velocities are not explicit in the dynamic and kinematic equations of rotation, the effect of the wind field on these equations is transitive through the aerodynamic moment coefficients. Even under wind perturbations, the equations are all the same as those without disturbance. Thus, it is not necessary to rewrite them again here.


### 2.1.3. Equations Derived Using Inertial Velocities under Wind Effect

All of the nonlinear equations under a wind condition derived using inertial velocities have the same form as those without perturbations, withV
→
e simply substituted for V
→. Thus, we can easily obtain the dynamic equations of force using the inertial velocities:(16)
u
˙
e
=
F
x
m
+
r
v
e
-
q
w
e
,
v
˙
e
=
F
y
m
+
p
w
e
-
r
u
e
,
w
˙
e
=
F
z
m
+
q
u
e
-
p
v
e
.Then, the equations of navigation are as follows:(17)
x
˙
e
=
u
e
cos
⁡
θ
cos
⁡
ψ
+
v
e
sin
⁡
ϕ
sin
⁡
θ
cos
⁡
ψ
-
cos
⁡
ϕ
sin
⁡
ψ
+
w
e
cos
⁡
ϕ
sin
⁡
θ
cos
⁡
ψ
+
sin
⁡
ϕ
sin
⁡
ψ
,
y
˙
e
=
u
e
cos
⁡
θ
sin
⁡
ψ
+
v
e
sin
⁡
ϕ
sin
⁡
θ
sin
⁡
ψ
+
cos
⁡
ϕ
cos
⁡
ψ
+
w
e
cos
⁡
ϕ
sin
⁡
θ
sin
⁡
ψ
-
sin
⁡
ϕ
cos
⁡
ψ
,
z
˙
e
=
-
u
e
sin
⁡
θ
+
v
e
sin
⁡
ψ
cos
⁡
θ
+
w
e
cos
⁡
ψ
cos
⁡
θ
.In accordance with the velocity triangle relationship, the components of the air-relative velocity vector in body coordinates are as follows:(18)
u
v
w
=
u
e
v
e
w
e
-
u
w
v
w
w
w
.Consequently, the airspeed and incidence can be obtained as follows:(19)
V
as
=
u
2
+
v
2
+
w
2
,
α
=
arctan
⁡
w
u
,
β
=
arcsin
⁡
v
V
as
.Differentiating (19), we can obtain the state equations of the flight parameters and find them the same as (14).


## 2.1.1. Nonlinear Equations without Perturbations

Typically, the forces acting on an aircraft primarily include the liftL, drag D, side force C, gravity G, and thrust T, which have the directions shown in Figure 1. Among these, the lift, drag, and side forces are called aerodynamic forces because they are caused by the motion of the aircraft relative to the air flow. Thus, their directions are defined in accordance with the direction of the air-relative velocity vector and generally computed in the air coordinates frame. In general, thrust is regarded as acting in the symmetric plane of the rigid aircraft with an inclination denoted by σ.Figure 1
Diagram of forces and velocities.LetV
→
=
u
v
w
T denote the air-relative velocity vector at the center of gravity, and let Ω
→
=
p
q
r
T be the angular velocity vector. Without perturbations, the inertial velocity is identical to the air-relative velocity. Then, the total force acting on the aircraft is represented in body coordinates as follows:(1)
F
→
=
F
x
F
y
F
z
=
T
cos
⁡
σ
0
T
sin
⁡
σ
+
C
B
A
-
D
C
-
L
+
C
B
E
0
0
m
g
,where the symbols C
B
A and C
B
E are the transform matrices from air coordinates and inertial coordinates to body coordinates, respectively. The minus sign in the drag and lift components means that their directions are the opposite to those of the coordinates adopted. According to Newton’s second law, the dynamic equation of force in body coordinates is as follows:(2)
m
V
→
˙
+
Ω
→
×
V
→
=
F
→
.Its components are rewritten as follows:(3)
u
˙
=
F
x
m
+
r
v
-
q
w
v
˙
=
F
y
m
+
p
w
-
r
u
w
˙
=
F
z
m
+
q
u
-
p
v
.LetV
as, α, β denote the airspeed, incidence, and sideslip angle, respectively, which correspond to the air-relative velocity by the following relationships:(4)
V
as
=
u
2
+
v
2
+
w
2
,
α
=
arctan
⁡
w
u
,
β
=
arcsin
⁡
v
V
as
.By differentiating the above with respect to time and replacingu
˙, v
˙, w
˙ with (3), we can obtain the equations of force in [14] as follows:(5)
V
˙
as
=
1
m
F
x
cos
⁡
α
cos
⁡
β
+
F
y
sin
⁡
β
+
F
z
sin
⁡
α
cos
⁡
β
,
α
˙
=
-
F
x
sin
⁡
α
+
F
z
cos
⁡
α
m
V
as
cos
⁡
β
+
q
-
tan
⁡
β
p
cos
⁡
α
+
r
sin
⁡
α
,
β
˙
=
1
m
V
as
-
F
x
cos
⁡
α
sin
⁡
β
+
F
y
cos
⁡
β
-
F
z
sin
⁡
α
sin
⁡
β
+
p
sin
⁡
α
-
r
cos
⁡
α
.Here, we can state that (5) is exactly the same as (3), in both theory and practice, as long as the forces are equal.Second, the equations of navigation are as follows:(6)
x
˙
=
u
cos
⁡
θ
cos
⁡
ψ
+
v
sin
⁡
ϕ
sin
⁡
θ
cos
⁡
ψ
-
cos
⁡
ϕ
sin
⁡
ψ
+
w
cos
⁡
ϕ
sin
⁡
θ
cos
⁡
ψ
+
sin
⁡
ϕ
sin
⁡
ψ
,
y
˙
=
u
cos
⁡
θ
sin
⁡
ψ
+
v
sin
⁡
ϕ
sin
⁡
θ
sin
⁡
ψ
+
cos
⁡
ϕ
cos
⁡
ψ
+
w
cos
⁡
ϕ
sin
⁡
θ
sin
⁡
ψ
-
sin
⁡
ϕ
cos
⁡
ψ
,
z
˙
=
-
u
sin
⁡
θ
+
v
sin
⁡
ψ
cos
⁡
θ
+
w
cos
⁡
ψ
cos
⁡
θ
.Third, the dynamic equations of moment in body coordinates are directly given without deduction:(7)
p
˙
=
1
I
x
x
I
z
z
-
I
x
z
2
I
x
x
-
I
y
y
+
I
z
z
I
x
z
p
q
+
I
y
y
I
z
z
-
I
z
z
2
-
I
x
z
2
q
r
+
I
z
z
l
+
I
x
z
N
,
q
˙
=
1
I
y
y
I
z
z
-
I
x
x
p
r
+
I
x
z
r
2
-
p
2
+
M
,
r
˙
=
1
I
x
x
I
z
z
-
I
x
z
2
I
x
x
2
-
I
x
x
I
y
y
+
I
x
z
2
p
q
+
-
I
x
x
+
I
y
y
-
I
z
z
I
x
z
q
r
+
I
x
z
l
+
I
x
x
N
.In the above,l, M, N are the components of the total external torque vector in body coordinates. I
x
x, I
y
y, I
z
z denote the moments of inertia, and I
x
z is the product of inertia of the rigid aircraft.Finally, the kinematic equations of rotation are as follows:(8)
ψ
˙
=
q
sin
⁡
ϕ
+
r
cos
⁡
ϕ
cos
⁡
θ
,
θ
˙
=
q
cos
⁡
ϕ
-
r
sin
⁡
ϕ
,
ϕ
˙
=
p
+
tan
⁡
θ
·
q
sin
⁡
ϕ
+
r
cos
⁡
ϕ
.


## 2.1.2. Equations Derived Using Air-Relative Velocities under Wind Effect

A vector triangle, namely,V
→
e
=
V
→
+
W
→, is the most fundamental relationship between velocities under wind. In order to distinguish these vectors, we use the subscript e to denote the vectors relative to the inertial coordinates and subscript w to denote the variables related to the wind velocity. Thus, let u
e, v
e, w
e denote the components of inertial velocity vector V
→
e in body coordinates, and let W
→
=
u
w
v
w
w
w
T denote the wind vector. To set up the nonlinear equations using air-relative velocities, we can substitute V
→
e in the relationship V
→
e
=
V
→
+
W
→ for the term V
→ in (2), as follows:(9)
m
V
→
˙
+
W
→
˙
+
Ω
→
×
V
→
+
W
→
=
F
→
.After arranging, we have the following:(10)
V
→
˙
=
F
→
m
-
Ω
→
×
V
→
-
W
→
˙
+
Ω
→
×
W
→
,where the term W
→
˙
+
Ω
→
×
W
→ is the absolute derivative of the wind velocity vector and called the effective force to be multiplied by the mass in [1, 12].The components of the above are as follows:(11)
u
˙
=
F
x
m
+
r
v
-
q
w
-
u
˙
w
+
r
v
w
-
q
w
w
,
v
˙
=
F
y
m
+
p
w
-
r
u
-
v
˙
w
+
p
w
w
-
r
u
w
,
w
˙
=
F
z
m
+
q
u
-
p
v
-
w
˙
w
+
q
u
w
-
p
v
w
.We can directly obtain the airspeed and incidence from the above through the mutual relationship in (4). The derivatives of the wind are as follows:(12)
u
˙
w
=
∂
u
w
∂
t
+
∂
u
w
∂
x
x
˙
e
+
∂
u
w
∂
y
y
˙
e
+
∂
u
w
∂
z
z
˙
e
,
v
˙
w
=
∂
v
w
∂
t
+
∂
v
w
∂
x
x
˙
e
+
∂
v
w
∂
y
y
˙
e
+
∂
v
w
∂
z
z
˙
e
,
w
˙
w
=
∂
w
w
∂
t
+
∂
w
w
∂
x
x
˙
e
+
∂
w
w
∂
y
y
˙
e
+
∂
w
w
∂
z
z
˙
e
.Since the motion of the air flow is slow enough relative to the airspeed of an aircraft, the wind velocity is generally regarded as steady with respect to time; namely,∂
u
w
/
∂
t
=
∂
v
w
/
∂
t
=
∂
w
w
/
∂
t
=
0. This is why we assume that the wind field satisfies the hypothesis of Taylor’s Frozen Field. The residual terms all include the wind gradients, as will be addressed later.Let(13)
F
→
w
=
F
w
x
F
w
y
F
w
z
=
m
u
˙
w
+
r
v
w
-
q
w
w
m
v
˙
w
+
p
w
w
-
r
u
w
m
w
˙
w
+
q
u
w
-
p
v
wdenote the effective force generated by the wind field. Then, the state equations for the airspeed, incidence, and sideslip angle in [14] are merely cited as follows:(14)
V
˙
as
=
1
m
F
x
-
F
w
x
cos
⁡
α
cos
⁡
β
+
F
y
-
F
w
y
sin
⁡
β
+
F
z
-
F
w
z
sin
⁡
α
cos
⁡
β
,
α
˙
=
-
F
x
-
F
w
x
sin
⁡
α
+
F
z
-
F
w
z
cos
⁡
α
m
V
as
cos
⁡
β
+
q
-
tan
⁡
β
p
cos
⁡
α
+
r
sin
⁡
α
,
β
˙
=
1
m
V
as
-
F
x
-
F
w
x
cos
⁡
α
sin
⁡
β
+
F
y
-
F
w
y
cos
⁡
β
-
F
z
-
F
w
z
sin
⁡
α
sin
⁡
β
+
p
sin
⁡
α
-
r
cos
⁡
α
.Compared with (5), the dynamic equations of force above are generally the same in form, only to subtract the effective force generated by the wind acceleration.Equations of navigation are obtained usingV
→
+
W
→ to replace V
→ in (6), as follows:(15)
x
˙
e
=
u
+
u
w
cos
⁡
θ
cos
⁡
ψ
+
v
+
v
w
sin
⁡
ϕ
sin
⁡
θ
cos
⁡
ψ
-
cos
⁡
ϕ
sin
⁡
ψ
+
w
+
w
w
cos
⁡
ϕ
sin
⁡
θ
cos
⁡
ψ
+
sin
⁡
ϕ
sin
⁡
ψ
,
y
˙
e
=
u
+
u
w
cos
⁡
θ
sin
⁡
ψ
+
v
+
v
w
sin
⁡
ϕ
sin
⁡
θ
sin
⁡
ψ
+
cos
⁡
ϕ
cos
⁡
ψ
+
w
+
w
w
cos
⁡
ϕ
sin
⁡
θ
sin
⁡
ψ
-
sin
⁡
ϕ
cos
⁡
ψ
,
z
˙
e
=
-
u
+
u
w
sin
⁡
θ
+
v
+
v
w
sin
⁡
ψ
cos
⁡
θ
+
w
+
w
w
cos
⁡
ψ
cos
⁡
θ
.Since the velocities are not explicit in the dynamic and kinematic equations of rotation, the effect of the wind field on these equations is transitive through the aerodynamic moment coefficients. Even under wind perturbations, the equations are all the same as those without disturbance. Thus, it is not necessary to rewrite them again here.


## 2.1.3. Equations Derived Using Inertial Velocities under Wind Effect

All of the nonlinear equations under a wind condition derived using inertial velocities have the same form as those without perturbations, withV
→
e simply substituted for V
→. Thus, we can easily obtain the dynamic equations of force using the inertial velocities:(16)
u
˙
e
=
F
x
m
+
r
v
e
-
q
w
e
,
v
˙
e
=
F
y
m
+
p
w
e
-
r
u
e
,
w
˙
e
=
F
z
m
+
q
u
e
-
p
v
e
.Then, the equations of navigation are as follows:(17)
x
˙
e
=
u
e
cos
⁡
θ
cos
⁡
ψ
+
v
e
sin
⁡
ϕ
sin
⁡
θ
cos
⁡
ψ
-
cos
⁡
ϕ
sin
⁡
ψ
+
w
e
cos
⁡
ϕ
sin
⁡
θ
cos
⁡
ψ
+
sin
⁡
ϕ
sin
⁡
ψ
,
y
˙
e
=
u
e
cos
⁡
θ
sin
⁡
ψ
+
v
e
sin
⁡
ϕ
sin
⁡
θ
sin
⁡
ψ
+
cos
⁡
ϕ
cos
⁡
ψ
+
w
e
cos
⁡
ϕ
sin
⁡
θ
sin
⁡
ψ
-
sin
⁡
ϕ
cos
⁡
ψ
,
z
˙
e
=
-
u
e
sin
⁡
θ
+
v
e
sin
⁡
ψ
cos
⁡
θ
+
w
e
cos
⁡
ψ
cos
⁡
θ
.In accordance with the velocity triangle relationship, the components of the air-relative velocity vector in body coordinates are as follows:(18)
u
v
w
=
u
e
v
e
w
e
-
u
w
v
w
w
w
.Consequently, the airspeed and incidence can be obtained as follows:(19)
V
as
=
u
2
+
v
2
+
w
2
,
α
=
arctan
⁡
w
u
,
β
=
arcsin
⁡
v
V
as
.Differentiating (19), we can obtain the state equations of the flight parameters and find them the same as (14).


## 2.2. Modeling of Aerodynamic Coefficients

Now, the data for modeling the aerodynamic coefficients can primarily be obtained from flight test records, wind tunnel experiment data, and aided computation methods. The nondimensional aerodynamic coefficients are typically dependent on such variables such as the incidence, sideslip angle, elevator deflection, rotational angular rates, rudder deflection, and aileron deflection. Of these, each aerodynamic derivative is given by a one- or two-dimensional table with respect to the current height and Mach and computed through linear interpolation when given the flight state parameters as inputs. The coefficients model adopted here uses components in body coordinates:(20)
C
x
=
C
x
0
+
C
x
α
α
+
C
x
δ
r
δ
r
,
C
y
=
C
y
β
β
+
C
y
α
β
α
β
+
C
y
δ
r
δ
r
+
C
y
δ
a
δ
a
,
C
z
=
C
z
α
α
+
C
z
δ
e
δ
e
,
C
l
=
C
l
β
β
+
C
l
β
α
α
β
+
C
l
p
p
b
2
V
+
C
l
r
r
b
2
V
+
C
l
δ
r
δ
r
+
C
l
δ
a
δ
a
+
C
l
δ
a
α
δ
a
α
,
C
m
=
C
m
0
+
C
m
α
α
+
C
m
q
q
c
-
V
+
C
m
δ
e
δ
e
+
C
m
r
r
b
2
V
,
C
n
=
C
n
0
+
C
n
β
α
α
β
+
C
n
β
β
+
C
n
p
p
b
2
V
+
C
n
r
r
b
2
V
+
C
n
δ
r
δ
r
+
C
n
δ
a
δ
a
.


## 2.3. Modeling of Wind Gradients

Since different parts of a rigid aircraft under a variable wind field may be influenced by different wind velocity vectors simultaneously, this dissimilarity would result in changes in the aerodynamic forces and moments acting on the aircraft [1]. Evans et al. [16] summarized the influence of spatial wind gradients upon the aerodynamic coefficients. In their thesis, a modified vortex-lattice method that took wind gradients into consideration was adopted to compute the aerodynamic coefficients of a wing and horizontal tail. Through their analyses and comparisons, they concluded that, under a microburst wind shear, wind gradients have such a significant influence upon the aerodynamic moments that they cannot be negligible. However, the impact on the aerodynamic forces is relatively unimportant. Actually, there are nine terms for wind gradients in spatial wind fields, and it is difficult to determine which is solely responsible for pitching, rolling, or yawing because of the interaction between these motions. In addition, the majority of them have no significant impact on the motion of an aircraft, especially under a small angle of attack. Etkin pointed out that a nonuniform distribution of air flow along the fuselage and wingspan would cause additional aerodynamic forces and moments. For the sake of simplification, based on the assumption that the wind velocity distribution is uniform along the coordinates, he gave four primary wind gradients, namely, ∂
w
w
/
∂
y, ∂
w
w
/
∂
x, ∂
v
w
/
∂
x, and ∂
u
w
/
∂
y. These will cause negative rolling, positive pitching, and negative and positive yawing, respectively, in body coordinates [17]. Thus, the rotation angular velocity used to compute the aerodynamic coefficients can be modified as follows:(21)
p
eff
q
eff
r
eff
=
p
q
r
+
-
∂
w
w
∂
y
∂
w
w
∂
x
∂
u
w
∂
y
-
∂
v
w
∂
x
.


## 3. Modeling of Microburst Wind Field

### 3.1. Vortex-Ring Principle

The principle of the vortex-ring model is briefly illuminated here based on [15]. In fact, the vertical component of the wind velocity must be zero on the ground boundary. To satisfy that, two vortex rings denoted by Γ and -
Γ, whose strengths are equal and opposite, can be arranged symmetrically about the ground. As shown in Figure 2, the ground coordinates O
x
y
z are fixed at an arbitrary position on the ground as its origin, and the center of the primary vortex ring is fixed at O
P
(
x
P
,
y
P
,
z
P
) in the ground coordinates. Similarly, the imaginary vortex ring is arranged at the symmetrical point O
I
(
x
P
,
y
P
,
-
z
P
). Any point N
(
x
N
,
y
N
,
z
N
) in the flow field is regarded as the reference point.Figure 2
Schematic diagram of vortex-ring model.LetR denote the radius of the vortex ring. According to the knowledge of fluid mechanics, the stream function is as follows:(22)
ψ
P
=
-
Γ
2
π
r
1
+
r
2
F
1
k
-
E
1
k
.In the above,Γ is the strength of the vortex pair; r
1 and r
2 denote the minimum and maximum distances from the reference point to the vortex ring, respectively; and F
1
(
k
) and E
1
(
k
) are the elliptic integral functions, where k
=
r
2
-
r
1
/
r
1
+
r
2 is generally called a scale factor. Practically, F
1
(
k
)
-
E
1
(
k
) can be replaced by its approximation: (23)
F
1
k
-
E
1
k
=
3.153
k
2
1
+
3
1
-
k
2
.In addition,r
1 and r
2 can be computed as follows:(24)
r
1
=
z
N
-
z
P
2
+
r
-
R
2
,
r
2
=
z
N
-
z
P
2
+
r
+
R
2
,where r is the radial distance to the central axis and computed as follows:(25)
r
=
x
N
-
x
P
2
+
y
N
-
y
P
2
.Then, we can obtain the stream function of the primary vortex ring:(26)
ψ
P
=
-
Γ
2
π
r
1
+
r
2
3.152
k
2
1
+
3
1
-
k
2
.In a similar way, the stream function of the imaginary vortex ring can be obtained:(27)
ψ
I
=
Γ
2
π
r
1
′
+
r
2
′
3.152
k
I
2
1
+
3
1
-
k
I
2
,where r
1
′, r
2
′, and k
I belong to the imaginary vortex ring, with the same meaning as r
1, r
2, and k, respectively.Sequentially, the stream function of the reference point can be obtained:(28)
ψ
=
ψ
P
+
ψ
I
.Finally, the induced velocity of this model is as follows:(29)
u
w
=
x
-
x
P
r
1
r
∂
ψ
∂
z
,
v
w
=
y
-
y
P
r
1
r
∂
ψ
∂
z
,
w
w
=
-
1
r
∂
ψ
∂
r
.From the above, we can easily find that the induced velocity is dependent on the distances relative to the vortex ring, radius, and strength of the vortex ring. In the meantime, we can also find the singularity that the velocity tends to infinity when the reference point approaches the filament of the vortex ring. In fact, there is a vortex core, due to viscosity, around the filament, and the velocity gradually decreases to zero in it. To solve this problem, a velocity damping factor was developed in [18]:(30)
ζ
=
1
-
exp
⁡
-
r
1
2
ε
d
2
,where r
1 is the distance in (24), ε is a weight coefficient used to adjust the degree of influence, and d is an empirical value related to the radius of the vortex core. The induced velocity can be computed using (29) directly outside the vortex core. In it, a velocity damping factor needs to be multiplied by the result. By means of these, the induced velocity will smoothly decrease to zero in the vortex core.


### 3.2. Double-Vortex-Ring Algorithm

The double-vortex-ring algorithm derived in [4] is used to build the microburst wind field here. On the basis of the vortex-ring model in Section 3.1, two vortex pairs with different parameters, including four vortex rings, are set symmetrical about the ground coordinates. The induced velocity and its damping factor for each vortex ring can be obtained by (29) and (30). Then, the whole velocity given in [4] is as follows:(31)
u
w
=
∏
i
=
1
4
ζ
i
∑
i
=
1
4
u
w
i
,
v
w
=
∏
i
=
1
4
ζ
i
∑
i
=
1
4
v
w
i
,
w
w
=
∏
i
=
1
4
ζ
i
∑
i
=
1
4
w
w
i
,where ∏
i
=
1
4
ζ
i, the product of the four damping factors, is called the total damping factor at the reference point.In [4], an algorithm was used to emulate the Dallas Fort Wind Shear encountered by two aircrafts, and its result was compared with the flight record data, finding it fit to simulate the microburst wind field. The parameters of the two vortex pairs, the large ring and small ring, are listed in Table 1. In computing the damping factors, ε is 0.1 and d is 2.7 times the radius of the vortex core.Table 1
Parameters of microburst wind field.
Parameter
Dimension
Large ring
Small ring
Ring radius
m
2134
396
Core radius
m
610
100
Vortex strength
m2/s
40131
12233
x position
m
760
90
y position
m
−90
1
z position
m
1036
245In addition, the profile of the microburst on the central axis is shown in Figure3. On the whole, the downdraft and head wind can apparently be seen.Figure 3
Profile of microburst on central axis.


## 3.1. Vortex-Ring Principle

The principle of the vortex-ring model is briefly illuminated here based on [15]. In fact, the vertical component of the wind velocity must be zero on the ground boundary. To satisfy that, two vortex rings denoted by Γ and -
Γ, whose strengths are equal and opposite, can be arranged symmetrically about the ground. As shown in Figure 2, the ground coordinates O
x
y
z are fixed at an arbitrary position on the ground as its origin, and the center of the primary vortex ring is fixed at O
P
(
x
P
,
y
P
,
z
P
) in the ground coordinates. Similarly, the imaginary vortex ring is arranged at the symmetrical point O
I
(
x
P
,
y
P
,
-
z
P
). Any point N
(
x
N
,
y
N
,
z
N
) in the flow field is regarded as the reference point.Figure 2
Schematic diagram of vortex-ring model.LetR denote the radius of the vortex ring. According to the knowledge of fluid mechanics, the stream function is as follows:(22)
ψ
P
=
-
Γ
2
π
r
1
+
r
2
F
1
k
-
E
1
k
.In the above,Γ is the strength of the vortex pair; r
1 and r
2 denote the minimum and maximum distances from the reference point to the vortex ring, respectively; and F
1
(
k
) and E
1
(
k
) are the elliptic integral functions, where k
=
r
2
-
r
1
/
r
1
+
r
2 is generally called a scale factor. Practically, F
1
(
k
)
-
E
1
(
k
) can be replaced by its approximation: (23)
F
1
k
-
E
1
k
=
3.153
k
2
1
+
3
1
-
k
2
.In addition,r
1 and r
2 can be computed as follows:(24)
r
1
=
z
N
-
z
P
2
+
r
-
R
2
,
r
2
=
z
N
-
z
P
2
+
r
+
R
2
,where r is the radial distance to the central axis and computed as follows:(25)
r
=
x
N
-
x
P
2
+
y
N
-
y
P
2
.Then, we can obtain the stream function of the primary vortex ring:(26)
ψ
P
=
-
Γ
2
π
r
1
+
r
2
3.152
k
2
1
+
3
1
-
k
2
.In a similar way, the stream function of the imaginary vortex ring can be obtained:(27)
ψ
I
=
Γ
2
π
r
1
′
+
r
2
′
3.152
k
I
2
1
+
3
1
-
k
I
2
,where r
1
′, r
2
′, and k
I belong to the imaginary vortex ring, with the same meaning as r
1, r
2, and k, respectively.Sequentially, the stream function of the reference point can be obtained:(28)
ψ
=
ψ
P
+
ψ
I
.Finally, the induced velocity of this model is as follows:(29)
u
w
=
x
-
x
P
r
1
r
∂
ψ
∂
z
,
v
w
=
y
-
y
P
r
1
r
∂
ψ
∂
z
,
w
w
=
-
1
r
∂
ψ
∂
r
.From the above, we can easily find that the induced velocity is dependent on the distances relative to the vortex ring, radius, and strength of the vortex ring. In the meantime, we can also find the singularity that the velocity tends to infinity when the reference point approaches the filament of the vortex ring. In fact, there is a vortex core, due to viscosity, around the filament, and the velocity gradually decreases to zero in it. To solve this problem, a velocity damping factor was developed in [18]:(30)
ζ
=
1
-
exp
⁡
-
r
1
2
ε
d
2
,where r
1 is the distance in (24), ε is a weight coefficient used to adjust the degree of influence, and d is an empirical value related to the radius of the vortex core. The induced velocity can be computed using (29) directly outside the vortex core. In it, a velocity damping factor needs to be multiplied by the result. By means of these, the induced velocity will smoothly decrease to zero in the vortex core.


## 3.2. Double-Vortex-Ring Algorithm

The double-vortex-ring algorithm derived in [4] is used to build the microburst wind field here. On the basis of the vortex-ring model in Section 3.1, two vortex pairs with different parameters, including four vortex rings, are set symmetrical about the ground coordinates. The induced velocity and its damping factor for each vortex ring can be obtained by (29) and (30). Then, the whole velocity given in [4] is as follows:(31)
u
w
=
∏
i
=
1
4
ζ
i
∑
i
=
1
4
u
w
i
,
v
w
=
∏
i
=
1
4
ζ
i
∑
i
=
1
4
v
w
i
,
w
w
=
∏
i
=
1
4
ζ
i
∑
i
=
1
4
w
w
i
,where ∏
i
=
1
4
ζ
i, the product of the four damping factors, is called the total damping factor at the reference point.In [4], an algorithm was used to emulate the Dallas Fort Wind Shear encountered by two aircrafts, and its result was compared with the flight record data, finding it fit to simulate the microburst wind field. The parameters of the two vortex pairs, the large ring and small ring, are listed in Table 1. In computing the damping factors, ε is 0.1 and d is 2.7 times the radius of the vortex core.Table 1
Parameters of microburst wind field.
Parameter
Dimension
Large ring
Small ring
Ring radius
m
2134
396
Core radius
m
610
100
Vortex strength
m2/s
40131
12233
x position
m
760
90
y position
m
−90
1
z position
m
1036
245In addition, the profile of the microburst on the central axis is shown in Figure3. On the whole, the downdraft and head wind can apparently be seen.Figure 3
Profile of microburst on central axis.


## 4. Simulation Experiments

Incorporating the microburst model of Section3, three simulation experiments are developed with a typical large aircraft that belongs to the C series of China. These three simulations utilize nonlinear equations without perturbations, dynamic equations derived using inertial velocities, and air-relative velocities. The configuration parameters of the aircraft are given in Table 2.Table 2
Parameters of aircraft.
Symbols
Value
Dimension
Illustration
m
7310
Kg
Aircraft mass
S
23.0
m2
Reference area
b
7.15
m
Wing span
c
4.002
m
Mean aerodynamic chord
I
x
x
450
Kg⋅m2
Inertia moment aroundx-axis
I
y
y
4400
Kg⋅m2
Inertia moment aroundy-axis
I
z
z
4700
Kg⋅m2
Inertia moment aroundz-axis
I
x
z
85.37
Kg⋅m2
Product of inertiaInitially, the aircraft is flying at 350 m with an airspeed of 132 m/s, heading north in the inertial frame without rolling and sideslip. Its incidence and pitch angle are equal to 9.17°, and its flight path angle is zero. The simulation platform will adopt a traditional fourth-order Runge-Kutta algorithm with a fixed step of 0.001 s.


### 4.1. Results of Simulation

In the simulation platform, the microburst wind velocity that the aircraft experiences is shown in Figure4. It produces a strong head wind of up to 20 m/s, which will undoubtedly have a significant influence on the motion of the aircraft.Figure 4
Wind velocity distribution.In Figure5, the blue lines denote the responses of the platform without perturbations, the black lines are the responses of the platform derived using inertial velocities, and the red lines are the responses of the platform using air-relative velocities. The responses of the four main variables, namely, the pitch angle, height, angle of attack, and airspeed, for the three platforms are compared. We can see that the magnitude of the responses of the platform using air-relative velocities fluctuates within a narrow range, while the responses of the platform using inertial velocities are relatively sharp. In addition, we get almost the same results using the simulink model in [14]. However, the changing trends of these are all the same like our results.Figure 5
Response of main variables.
(a)
(b)
(c)
(d)The changes in motion must be caused by the total forces and moments. Thus, in Figure6, the effective force and total force of the air-relative model are shown and compared. The initiation flight states are almost trimmed, so the external forces acting on the aircraft approach zero. Thus, the effective force is nearly equal to the total force. As has been shown in Section 2, (14) can be derived using both air-relative and inertial velocities, and thus we can deduce that the two forms of nonlinear equations of motion are identical in nature. Since the differences in motion lie in the total force acting on the rigid body, at the root, the estimation of the effective force induced by the wind velocity is too simple to approximate the changes in the aerodynamic forces and moments. Considering it multiplied by a weight factor, its changing trends will be maintained, and the responses of the two platforms will stand much closer, as shown in Figure 7.Figure 6
Effective force and total force.Figure 7
Responses after effective force is modified.
(a)
(b)
(c)
(d)


### 4.2. Discussions and Analyses

(1)
In these numerical simulation experiments, through the responses of the four main state variables of the nonlinear equations derived using the air-relative velocities and inertial velocities, the differences between them are primarily compared and analyzed. When the wind velocity is entered, the equations derived using the inertial velocities promptly update the airspeed using (18) and (19) and then compute all the aerodynamic coefficients using the changed airspeed. The equations derived using the air-relative velocities first estimate the effective force generated by the wind velocity and then compute the aerodynamic coefficients using the original airspeed. Obviously, the relationship between the aerodynamic coefficients and state variables (such as the airspeed and incidence) is seriously nonlinear. It is impossible to use the estimation of the effective force to exactly compensate for the different values between the aerodynamic forces computed by the changed airspeed and the original one. Once the inputs change, all of the state variables will change because of the interaction between all the motions of the aircraft.
(2)
Given that the aerodynamic coefficients are defined in accordance with the air coordinates, in theory, there is no doubt that we generally prefer using air-relative velocities to compute and analyze aerodynamic forces. However, in realistic situations, once the wind velocity vector is entered, the attitude of the aircraft is held, and the incidence is changed at the same time. This process is exactly reflected using the inertial velocities. In the platform using the effective forces constructed in [14], the responses of the state variables, without perturbations and under an intense wind field, are nearly the same. This result violates the characteristics of motion under wind, which in turn demonstrates the deviation using the air-relative velocities. Anyway, these results are not compared with realistic airplane because the data is not yet acquired. As such, we expect to indicate the percentage error in the upcoming works.


## 4.1. Results of Simulation

In the simulation platform, the microburst wind velocity that the aircraft experiences is shown in Figure4. It produces a strong head wind of up to 20 m/s, which will undoubtedly have a significant influence on the motion of the aircraft.Figure 4
Wind velocity distribution.In Figure5, the blue lines denote the responses of the platform without perturbations, the black lines are the responses of the platform derived using inertial velocities, and the red lines are the responses of the platform using air-relative velocities. The responses of the four main variables, namely, the pitch angle, height, angle of attack, and airspeed, for the three platforms are compared. We can see that the magnitude of the responses of the platform using air-relative velocities fluctuates within a narrow range, while the responses of the platform using inertial velocities are relatively sharp. In addition, we get almost the same results using the simulink model in [14]. However, the changing trends of these are all the same like our results.Figure 5
Response of main variables.
(a)
(b)
(c)
(d)The changes in motion must be caused by the total forces and moments. Thus, in Figure6, the effective force and total force of the air-relative model are shown and compared. The initiation flight states are almost trimmed, so the external forces acting on the aircraft approach zero. Thus, the effective force is nearly equal to the total force. As has been shown in Section 2, (14) can be derived using both air-relative and inertial velocities, and thus we can deduce that the two forms of nonlinear equations of motion are identical in nature. Since the differences in motion lie in the total force acting on the rigid body, at the root, the estimation of the effective force induced by the wind velocity is too simple to approximate the changes in the aerodynamic forces and moments. Considering it multiplied by a weight factor, its changing trends will be maintained, and the responses of the two platforms will stand much closer, as shown in Figure 7.Figure 6
Effective force and total force.Figure 7
Responses after effective force is modified.
(a)
(b)
(c)
(d)


## 4.2. Discussions and Analyses

(1)
In these numerical simulation experiments, through the responses of the four main state variables of the nonlinear equations derived using the air-relative velocities and inertial velocities, the differences between them are primarily compared and analyzed. When the wind velocity is entered, the equations derived using the inertial velocities promptly update the airspeed using (18) and (19) and then compute all the aerodynamic coefficients using the changed airspeed. The equations derived using the air-relative velocities first estimate the effective force generated by the wind velocity and then compute the aerodynamic coefficients using the original airspeed. Obviously, the relationship between the aerodynamic coefficients and state variables (such as the airspeed and incidence) is seriously nonlinear. It is impossible to use the estimation of the effective force to exactly compensate for the different values between the aerodynamic forces computed by the changed airspeed and the original one. Once the inputs change, all of the state variables will change because of the interaction between all the motions of the aircraft.
(2)
Given that the aerodynamic coefficients are defined in accordance with the air coordinates, in theory, there is no doubt that we generally prefer using air-relative velocities to compute and analyze aerodynamic forces. However, in realistic situations, once the wind velocity vector is entered, the attitude of the aircraft is held, and the incidence is changed at the same time. This process is exactly reflected using the inertial velocities. In the platform using the effective forces constructed in [14], the responses of the state variables, without perturbations and under an intense wind field, are nearly the same. This result violates the characteristics of motion under wind, which in turn demonstrates the deviation using the air-relative velocities. Anyway, these results are not compared with realistic airplane because the data is not yet acquired. As such, we expect to indicate the percentage error in the upcoming works.


## 5. Conclusions

This paper analyzed and compared the nonlinear equations of a rigid aircraft derived using the air-relative velocities and inertial velocities under a wind effect. By deducing the dynamic equations for the airspeed, incidence, and sideslip angle, we found them to be identical in nature. In addition, the results of the simulation experiments showed that the tendencies of responses from the two models were almost consistent. However, the response amplitude when using the air-relative velocities was inconspicuous even under a fierce wind field, which was a result of the estimation error of the effective wind force. Thus, this paper showed a preference for the equations derived using inertial velocities because of their convenience and accuracy in practice. We will make further efforts to improve methods to estimate the effective wind force, which we believe will provide an intuitive explanation for the wind impact on aircrafts.



---

*Source: 2371594-2016-06-30.xml*