# Coordinated Control of Pressure Difference and Rising Velocity for Stratospheric Airship with Thermal Effects

Authors: Luhe Hong; Hui-Yu Jin; Xianwu Lin; Weiyao Lan
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Ascending control of stratospheric airship is a challenging control problem, especially if both the rising velocity and
the pressure difference between the inside and outside of the airship are required to be controlled simultaneously during ascending.
In this paper, a coordinated scheme to control pressure difference and rising velocity of stratospheric airship with vector thrust
is presented. With the control scheme, the airship maintains the pressure difference by exhausting air with feedback control. At
the same time, the supplemental thrust is generated to compensate the buoyancy fluctuation caused by exhausting air so that the
airship’s vertical velocity can track a given reference trajectory. Simulations show that the coordinated control scheme ensures that
the airship rises to the altitude of 20 km steadily and rapidly while the pressure difference is always in the safe range. Furthermore,
the control scheme is robust enough to the thermal disturbance caused by solar radiation and other thermal processes, which is
calculated with partial differential equations.


---

## Body

## 1. Introduction

Stratospheric airship, a low-speed near space aerocraft, mainly relies on static buoyancy to ascend and maintain itself in working height [1]. Because of its unparalleled advantages such as low power consumption, high security, and large shipping capacity, the prospect of airship has drawn increasing concern around the world, while airship technology has become a cutting-edge and hot issue for academic research [2].Rising to the resident height safely is a prerequisite for stratospheric airship to work properly [3]. During ascending, the stratospheric airship exhausts inner air continuously to guarantee that the pressure inside the airship is slightly greater than the atmospheric pressure outside. If the internal pressure is not enough, the airship balloon may be deflated by external atmospheric pressure, resulting in an overall structural deformation. If the internal pressure is too large, the balloon may be in excessive tension even rupture. Therefore, to control the pressure difference between inside and outside of the airship is extremely important throughout the whole process of ascending [4].From the perspective of control theory, pressure difference control during ascending is a difficult problem. It requires simultaneously high accuracy, rapidity, and robustness. First, the pressure difference of airship must be guaranteed in a very limited safe range. For example, a typical value of the internal pressure greater than the external pressure is about 300 Pa~600 Pa [5], which means the allowable maximum fluctuating value is only 300 Pa and less than 3‰ of the internal pressure of airship on the ground. To guarantee the pressure difference is always in the safe range, the control approach must be highly accurate. Second, in practice, the airship must have an appropriate increasing velocity, which leads to quickly decreasing the atmosphere pressure outside of the airship. For instance, if the airship has altitude 1 km and rising velocity 1 m/s, its external atmosphere pressure decreases by approximately 11 Pa per second. To maintain pressure difference in the safe range, the internal pressure must be adjusted quickly enough. Third and most important, in rising process there exist nonignorable unmodeled dynamics. Thermal effects, such as solar radiation and force convection, can significantly affect the internal pressure of the airship and the pressure difference. Unfortunately, these thermal effects are too complicated to be modeled with precise and brief mathematical models used for online control. So, the pressure difference control has to view these thermal effects as unmodeled dynamics and to be robust enough.In the current literature, the pressure difference control during ascending has not been carefully addressed. Seldom research considers pressure difference, rising performances, and thermal effects simultaneously. Affected by low-altitude airship, many researchers focus only on the flight dynamic performance, such as Zheng et al. [6], Bestaoui and Kahale [7], Zhang et al. [8], and Mueller et al. [9]. In [10, 11], two approaches are used, respectively, to control the pressure difference during ascending while the rising performances such as velocity and smoothness are not guaranteed. A coordinated control scheme of pressure difference and height of stratospheric airship was proposed by Wu et al. [4], but with the scheme the airship may oscillate in the vertical direction. On the other hand, Guo and Zhu [12] and Shi et al. [13] considered the thermodynamic interference during airship ascending, while the pressure difference and flight dynamics were not considered.This paper presents a new ascending control scheme for stratospheric airship which has vector thrust. By exhausting air and the vector thrust, the pressure difference and rising velocity of the airship are coordinately controlled so that the airship can rise to the resident height steadily and rapidly. In the scheme, the priority is given to control pressure difference by exhausting air with feedback. Then a controller is designed by pole placement method to adjust the thrust to compensate the rising velocity change caused by exhausting air and to track a given altitude-velocity reference. A thermal model of the airship is established to estimate thermal effects caused by direct solar radiation and other disturbances and to verify the robustness of the control scheme. Simulations show that, with the coordinated control, the airship can rise to the altitude of 20 km within 8 hours against thermal disturbance. During ascending, the pressure difference is always controlled in the range of 340 Pa~420 Pa while the airship does not oscillate in the vertical direction.The remainder of the paper is organized as follows. Section2 presents the airship features. In Section 3, the mathematical models of environment and airship are established. In Section 4, the corresponding coordinated control strategy is given. Section 5 gives the simulation results and analysis while Section 6 is for concluding remarks.


## 2. Airship Features

According to the different structures, traditional airships can be divided into three types: blimp airship, semirigid airship, and rigid airship. This paper investigates semirigid airship. A rigid dragon skeleton, which can maintain good structural shape and withstand a certain twisted moment, is designed in the lower part of the airship. Airship hull is designed as water droplets form, synthesized based on ellipsoid approximation, and it is the same minor axis length between two ellipsoids. The lengths of the short axis and the two long axes are denoted asb, a
1, and a
2, respectively. The total length of airship is L. Inside airship, there are helium bags and air ballonets. Some exhaust valves and inspiratory valves controlled by onboard computer are equipped in the air ballonets. Thrust devices equipped in the abdominal and tail of airship are also controlled by onboard computer. The overall structure schematic diagram of airship is shown in Figure 1.Figure 1
Structure schematic diagram of airship.Under normal conditions, helium bags are sealed and they are constant without mass exchanging with the outside. Helium bags transfer pressure to air ballonets through diaphragm, so the pressure remains the same everywhere in the internal of airship. To maintain the airship geometry and ensure the pressure within the safe range, the air ballonets exchange mass with outside air through valve or blower to maintain a certain pressure difference between the inside and outside of airship. Therefore, helium bags and air ballonets of stratospheric airship are divided into two control bodies for analysis.The total mass of airship is(1)
M
=
m
a
i
r
+
m
H
e
+
m
s
,where m
a
i
r is the mass of air, m
H
e is the mass of helium, and m
s is the mass of payload and ballonet material. m
H
e and m
s are constant, so the total mass takes into account m
a
i
r that is necessarily dragged along with the balloon.The volume of airship is(2)
V
=
2
3
π
a
1
+
a
2
b
2
,
V
=
V
a
i
r
+
V
H
e
,where V
a
i
r and V
H
e are the volume of internal air and helium, respectively.Because the pressure in the internal of airship is equal, it can be expressed as(3)
P
i
n
=
P
a
i
r
=
P
H
e
,where P
i
n is the internal pressure of airship and P
a
i
r and P
H
e are pressure of internal air and helium, respectively. The pressure difference between the internal and external of airship is(4)
Δ
P
=
P
i
n
-
P
a
t
m
,where P
a
t
m means the atmosphere pressure.Let the take-off point of airship beh
0
=
0 and let the final working height be h
f
=
20
k
m. When the airship hovers itself in the working height for a long time, the inner temperature of gas can be assumed to be a constant which equals the atmosphere temperature at h
f
=
20
k
m. The inner gas obeys Boyle’s Law; that is,(5)
ρ
H
e
h
f
=
P
a
t
m
h
f
+
Δ
P
R
H
e
T
a
t
m
h
f
,
ρ
a
i
r
h
f
=
P
a
t
m
h
f
+
Δ
P
R
a
i
r
T
a
t
m
h
f
,where T
a
t
m is the temperature of the atmosphere, ρ
a
i
r and ρ
H
e are the density of internal air and helium, respectively, and R
a
i
r and R
H
e are the gas constant of air and helium, respectively. The final volume of helium is(6)
V
H
e
h
f
=
m
H
e
ρ
H
e
h
f
.Therefore, the final mass of internal air is(7)
m
a
i
r
h
f
=
ρ
a
i
r
h
f
V
-
m
H
e
ρ
H
e
h
f
.Buoyancy force F
b
=
ρ
a
t
m
g
V is equal to the weight of airship G
=
M
g at h
f, so we have(8)
F
b
h
f
=
G
h
f
,
ρ
a
t
m
h
f
g
V
=
M
h
f
g
=
m
a
i
r
h
f
+
m
H
e
+
m
s
g
,where g is the gravitational acceleration. Thereby, m
s is(9)
m
s
=
ρ
a
t
m
h
f
V
-
ρ
a
i
r
h
f
V
-
m
H
e
ρ
H
e
h
f
-
m
H
e
.Because m
s is fixed and known, m
H
e can be solved by (9). Therefore, the initial mass of internal air is(10)
m
a
i
r
h
0
=
ρ
a
i
r
h
0
V
-
m
H
e
ρ
H
e
h
0
.


## 3. Mathematical Models of Environment and Airship

Coordinated control of pressure difference and rising velocity for stratospheric airship is closely related to atmospheric environment, solar radiation, and the motion state of airship during ascending, so it is necessary to establish the relevant models.


### 3.1. Atmosphere Model

Based on the standard atmosphere, the ambient air density and temperature are modeled as functions of the altitude. They can be expressed as follows [14]:(11)
ρ
a
t
m
h
=
ρ
a
t
m
0
T
a
t
m
h
T
a
t
m
0
n
1
,
h
0
≤
h
≤
h
m
;
ρ
a
t
m
h
m
e
n
2
h
-
h
m
,
h
m
<
h
≤
h
f
,
T
a
t
m
h
=
T
a
t
m
0
+
d
T
·
h
,
h
0
≤
h
≤
h
m
;
216.65
K
,
h
m
<
h
≤
h
f
,where n
1
=
4.256, n
2
=
-
1.578
×
1
0
-
4, h
m
=
11000
m, ρ
a
t
m
(
0
)
=
1.225
k
g
/
m
3, ρ
a
t
m
(
11000
)
=
0.3639
k
g
/
m
3, T
a
t
m
(
0
)
=
288.15
K, and d
T
=
-
0.0065
K.By ideal gas law, the ambient air satisfies(12)
P
a
t
m
h
=
ρ
a
t
m
h
R
a
i
r
T
a
t
m
h
.According to (11), there is(13)
d
ρ
a
t
m
h
d
h
=
-
a
0
ρ
a
t
m
h
,where(14)
a
0
=
0.96
×
1
0
-
4
,
h
0
≤
h
≤
h
m
;
1.578
×
1
0
-
4
,
h
m
<
h
≤
h
f
.Combining (11)~(13), the change rate of atmosphere pressure P
˙
a
t
m can be written as(15)
P
˙
a
t
m
=
ρ
a
t
m
h
R
a
i
r
v
d
T
-
a
0
T
a
t
m
h
,
h
0
≤
h
≤
h
m
;
-
a
0
ρ
a
t
m
h
R
a
i
r
v
T
a
t
m
h
,
h
m
<
h
≤
h
f
,where v is the vertical velocity of airship. Therefore, P
˙
a
t
m is related to the flight altitude h and velocity v of airship.


### 3.2. Solar Model

The methods calculating the value of solar radiation are studied by Ran et al. [15], Farley [16], and Dai et al. [14]. The direct solar irradiance flux I
D and the diffuse irradiance flux I
S can be expressed as(16)
I
D
=
I
0
1
+
e
cos
⁡
ζ
1
-
e
2
2
p
t
m
,
I
S
=
1
2
I
D
sin
⁡
ω
1
-
p
t
m
1
-
1.4
ln
⁡
p
t
,where I
0
=
1367
W
/
m
2 is solar constant, e is the orbital eccentricity, for earth e
=
1.671
×
1
0
-
3, and the value of 0.6~0.7 is adopted for atmospheric transmittance p
t. The air mass ratio m and the true anomaly ζ can be calculated as follows:(17)
m
=
P
a
t
m
P
0
1229
+
614
sin
⁡
ω
2
-
614
sin
⁡
ω
,
ζ
=
M
A
+
2
e
sin
⁡
M
A
+
1.25
e
2
sin
⁡
2
M
A
,where P
0 is the air pressure on the sea level, ω is the solar elevation angle (
ω
≥
0
), and M
A
=
2
n
π
/
365, where n is the day number in a year. For example, n
=
1 denotes January 1.The reflected solar irradiance fluxI
R also increases the irradiance flux, which is given as(18)
I
R
=
C
R
I
D
sin
⁡
ω
+
I
S
,where C
R is the albedo factor.The total solar irradiance flux can be expressed as(19)
I
=
I
D
+
I
S
+
I
R
.


### 3.3. Kinematic and Dynamic Models of Airship

The airship’s motion is subject to the net buoyancyB of airship and the atmosphere wind v
w. It can be assumed that the horizontal velocity of airship is equal to the atmospheric wind [17]. Therefore, the horizontal motion of airship can be assumed at the wind speed. An analysis of the vertical force balance on the airship is shown in Figure 2.Figure 2
Analysis of the vertical force balance on the airship.The thrust of airshipT is generated by vector thrust mechanism.The net buoyancy of airship is(20)
B
=
F
b
-
G
.The drag force of airshipD can be calculated by(21)
D
=
1
2
ρ
a
t
m
h
v
v
S
C
d
,where S
=
V
2
/
3 is the reference area and C
d is the drag coefficient.For the range of Reynolds numbers Re up to 107, the relationship between the drag coefficient C
d and Reynolds numbers Re was reported by Almedeij [18]. When Reynolds numbers are greater than 106, the drag coefficient is stable at the value of about 0.1, while for Reynolds numbers lower than 105 the constant drag coefficient of 0.5 can be used. A drag coefficient model for the ascent balloon in a Reynolds number range of 105~106 is explored by Conner and Arena [19], which can be expressed as follows:(22)
C
d
=
0.72
-
2.57
×
1
0
-
6
R
e
+
4.71
×
1
0
-
12
R
e
2
-
4.04
×
1
0
-
18
R
e
3
+
1.31
×
1
0
-
24
R
e
.The basic equations of airship during ascending can be expressed as follows:(23)
h
˙
=
v
,
v
˙
=
B
-
D
+
T
M
,
M
˙
=
m
˙
a
i
r
,where m
˙
a
i
r is the change rate of internal air mass. It is assumed that the airship only exhausts air and does not inhale air in this paper, so m
˙
a
i
r
≤
0.


### 3.4. Thermal Effects of Airship

The thermal effects of airship come from external and internal environment. External environment includes direct solar radiation, diffuse solar radiation, reflected radiation, infrared radiation, and forced convection between atmosphere and film. Internal environment includes infrared radiation among inner surface area elements and natural convection between inner gas and film. Thermal environment of airship is shown in Figure3.Figure 3
Thermal environment of airship.


#### 3.4.1. Membrane Temperature Model

The surface of airship is divided intoM
×
N grids. The transient energy equation for each grid A
i
j can be expressed as follows [20]:(24)
d
s
ρ
s
c
T
˙
i
j
=
k
a
t
m
,
i
j
T
a
t
m
-
T
i
j
+
k
i
n
,
i
j
T
i
n
-
T
i
j
-
ϵ
1
+
ϵ
2
σ
T
i
j
4
+
∑
k
=
1
4
q
k
,
i
j
+
∑
I
=
1
M
∑
J
=
1
N
ϵ
2
E
I
J
A
i
j
·
∫
A
i
j
∫
A
I
J
cos
⁡
α
i
j
cos
⁡
α
I
J
π
l
2
d
A
i
j
d
A
I
J
,where d
s is the thickness of the film, ρ and c are the density and specific heat of the film, respectively, T
i
j is the temperature of element A
i
j, T
i
n is the inner gas temperature which contacts to the film, k
a
t
m
,
i
j is the heat transfer coefficient between atmosphere and film, k
i
n
,
i
j is the heat transfer coefficient between inner gas and film, ϵ
1 and ϵ
2 are the emissivity of the internal and external of film, respectively, E
I
J is the radiation heat flux of element A
I
J, α
i
j and α
I
J are the included angle between the internal normal vector n
→
i
j, n
→
I
J, and the connection of two elements, respectively, and q
k
,
i
j is the absorbed radiation heat flux of external film, including direct solar radiation q
1
,
i
j, diffuse solar radiation q
2
,
i
j, reflected radiation q
3
,
i
j, and infrared radiation q
4
,
i
j. The absorbed solar radiation [21] and absorbed infrared radiation [13] can be expressed as follows:(i)
The absorbed direct solar radiationq
1
,
i
j can be expressed as(25)
q
1
,
i
j
=
α
s
A
i
j
I
D
cos
⁡
β
i
j
,
whereα
s is the solar absorptivity of the film and β
i
j is the included angle between the element normal and the solar irradiation.
(ii)
The absorbed diffuse solar radiationq
2
,
i
j can be expressed as(26)
q
2
,
i
j
=
α
s
A
i
j
I
S
0.5
-
0.5
cos
⁡
ϕ
i
j
,
whereϕ
i
j is the included angle between the element normal and the gravity direction.
(iii)
The absorbed reflected solar radiationq
3
,
i
j can be expressed as(27)
q
3
,
i
j
=
α
s
A
i
j
I
R
0.5
+
0.5
cos
⁡
ϕ
i
j
.
(iv)
The absorbed infrared radiationq
4
,
i
j can be expressed as(28)
q
4
,
i
j
=
ϵ
2
A
i
j
σ
T
a
t
m
4
-
T
i
j
4
,
whereσ is the Stefan-Boltzmann constant.


#### 3.4.2. Energy Equations of Inner Gas

Supposing that the temperature in one balloon is the same, according to the First Law of Thermodynamics, there are(29)
C
p
,
a
i
r
m
a
i
r
T
˙
a
i
r
=
Q
a
i
r
+
V
a
i
r
P
˙
i
n
-
C
p
,
a
i
r
T
a
i
r
m
˙
a
i
r
,
C
p
,
H
e
m
H
e
T
˙
H
e
=
Q
H
e
+
V
H
e
P
˙
i
n
,where C
p
,
a
i
r and C
p
,
H
e are the isobaric specific heats of air and helium, respectively. T
a
i
r and T
H
e are the temperature of air and helium, respectively. Q
a
i
r and Q
H
e are the quantity of heat convection of air and helium, respectively, which are given by(30)
Q
a
i
r
=
∑
i
=
1
M
∑
j
=
1
N
k
a
i
r
,
i
j
T
i
j
-
T
a
i
r
A
i
j
+
k
a
i
r
,
d
i
a
T
d
i
a
-
T
a
i
r
S
d
i
a
,
a
i
r
,
Q
H
e
=
∑
i
=
1
M
∑
j
=
1
N
k
H
e
,
i
j
T
i
j
-
T
H
e
A
i
j
+
k
H
e
,
d
i
a
T
d
i
a
-
T
H
e
S
d
i
a
,
H
e
,where T
d
i
a is the temperature of diaphragm between helium bag and air ballonet. k
a
i
r
,
i
j and k
a
i
r
,
d
i
a are the convection heat transfer coefficient between and film, diaphragm, and internal air. k
a
i
r
,
i
j, k
a
i
r
,
d
i
a, k
H
e
,
i
j and k
H
e
,
d
i
a are the convection heat transfer coefficients between internal air and film, internal air and diaphragm, helium and film, helium and diaphragm respectively. S
a
i
r
,
d
i
a and S
H
e
,
d
i
a are the contact areas between internal air and diaphragm, helium and diaphragm respectively. P
˙
i
n is the change rate of internal pressure which is given by(31)
V
P
˙
i
n
=
m
a
i
r
R
a
i
r
T
˙
a
i
r
+
m
H
e
R
H
e
T
˙
H
e
-
m
˙
a
i
r
R
a
i
r
T
a
i
r
.Combining (29)~(31), the change rate of inner air, helium, and internal pressure can be expressed as(32)
T
˙
a
i
r
=
Q
a
i
r
R
H
e
V
H
e
-
Q
H
e
R
H
e
V
a
i
r
-
C
p
,
H
e
Q
a
i
r
V
+
C
p
,
H
e
R
a
i
r
T
a
i
r
V
a
i
r
m
˙
a
i
r
-
C
p
,
a
i
r
R
H
e
T
a
i
r
V
H
e
m
˙
a
i
r
+
C
p
,
H
e
C
p
,
a
i
r
T
a
i
r
V
m
˙
a
i
r
m
a
i
r
C
p
,
H
e
R
a
i
r
V
a
i
r
+
C
p
,
a
i
r
R
H
e
V
H
e
-
C
p
,
H
e
C
p
,
a
i
r
V
,
T
˙
H
e
=
Q
H
e
R
a
i
r
V
a
i
r
-
Q
a
i
r
R
a
i
r
V
H
e
-
C
p
,
a
i
r
Q
H
e
V
+
2
C
p
,
a
i
r
R
a
i
r
T
a
i
r
V
H
e
m
˙
a
i
r
m
H
e
C
p
,
H
e
R
a
i
r
V
a
i
r
+
C
p
,
a
i
r
R
H
e
V
H
e
-
C
p
,
H
e
C
p
,
a
i
r
V
,
P
˙
i
n
=
-
C
p
,
H
e
Q
a
i
r
R
a
i
r
+
C
p
,
a
i
r
Q
H
e
R
H
e
-
2
C
p
,
H
e
C
p
,
a
i
r
R
a
i
r
T
a
i
r
m
˙
a
i
r
C
p
,
H
e
R
a
i
r
V
a
i
r
+
C
p
,
a
i
r
R
H
e
V
H
e
-
C
p
,
H
e
C
p
,
a
i
r
V
.


### 3.5. Actuator Model

Actuator includes exhaust valve and vector thrust mechanism. The exhaust valve expels inner air to maintain the pressure difference while vector thrust mechanism generates the vertical thrust to help control rising velocity.The parameters of exhaust valve are as follows: the diameter of valve isd
a
i
r, the adjustable opening area is S
a
i
r
=
π
d
a
i
r, the opening degree is f
a
i
r
(
t
), and flow resistance is C
a
i
r. According to Bernoulli’s principle, there is(33)
m
˙
a
i
r
t
=
2
P
i
n
t
Δ
P
t
R
a
i
r
T
a
i
r
t
S
a
i
r
C
a
i
r
f
a
i
r
t
.The effect of exhausting air is relevant not only to opening degree but also to the current thermodynamic states of inner gas.The vector thrust mechanism can continuously generate the vertical thrustT, and the range of T is 0~1000 N.


## 3.1. Atmosphere Model

Based on the standard atmosphere, the ambient air density and temperature are modeled as functions of the altitude. They can be expressed as follows [14]:(11)
ρ
a
t
m
h
=
ρ
a
t
m
0
T
a
t
m
h
T
a
t
m
0
n
1
,
h
0
≤
h
≤
h
m
;
ρ
a
t
m
h
m
e
n
2
h
-
h
m
,
h
m
<
h
≤
h
f
,
T
a
t
m
h
=
T
a
t
m
0
+
d
T
·
h
,
h
0
≤
h
≤
h
m
;
216.65
K
,
h
m
<
h
≤
h
f
,where n
1
=
4.256, n
2
=
-
1.578
×
1
0
-
4, h
m
=
11000
m, ρ
a
t
m
(
0
)
=
1.225
k
g
/
m
3, ρ
a
t
m
(
11000
)
=
0.3639
k
g
/
m
3, T
a
t
m
(
0
)
=
288.15
K, and d
T
=
-
0.0065
K.By ideal gas law, the ambient air satisfies(12)
P
a
t
m
h
=
ρ
a
t
m
h
R
a
i
r
T
a
t
m
h
.According to (11), there is(13)
d
ρ
a
t
m
h
d
h
=
-
a
0
ρ
a
t
m
h
,where(14)
a
0
=
0.96
×
1
0
-
4
,
h
0
≤
h
≤
h
m
;
1.578
×
1
0
-
4
,
h
m
<
h
≤
h
f
.Combining (11)~(13), the change rate of atmosphere pressure P
˙
a
t
m can be written as(15)
P
˙
a
t
m
=
ρ
a
t
m
h
R
a
i
r
v
d
T
-
a
0
T
a
t
m
h
,
h
0
≤
h
≤
h
m
;
-
a
0
ρ
a
t
m
h
R
a
i
r
v
T
a
t
m
h
,
h
m
<
h
≤
h
f
,where v is the vertical velocity of airship. Therefore, P
˙
a
t
m is related to the flight altitude h and velocity v of airship.


## 3.2. Solar Model

The methods calculating the value of solar radiation are studied by Ran et al. [15], Farley [16], and Dai et al. [14]. The direct solar irradiance flux I
D and the diffuse irradiance flux I
S can be expressed as(16)
I
D
=
I
0
1
+
e
cos
⁡
ζ
1
-
e
2
2
p
t
m
,
I
S
=
1
2
I
D
sin
⁡
ω
1
-
p
t
m
1
-
1.4
ln
⁡
p
t
,where I
0
=
1367
W
/
m
2 is solar constant, e is the orbital eccentricity, for earth e
=
1.671
×
1
0
-
3, and the value of 0.6~0.7 is adopted for atmospheric transmittance p
t. The air mass ratio m and the true anomaly ζ can be calculated as follows:(17)
m
=
P
a
t
m
P
0
1229
+
614
sin
⁡
ω
2
-
614
sin
⁡
ω
,
ζ
=
M
A
+
2
e
sin
⁡
M
A
+
1.25
e
2
sin
⁡
2
M
A
,where P
0 is the air pressure on the sea level, ω is the solar elevation angle (
ω
≥
0
), and M
A
=
2
n
π
/
365, where n is the day number in a year. For example, n
=
1 denotes January 1.The reflected solar irradiance fluxI
R also increases the irradiance flux, which is given as(18)
I
R
=
C
R
I
D
sin
⁡
ω
+
I
S
,where C
R is the albedo factor.The total solar irradiance flux can be expressed as(19)
I
=
I
D
+
I
S
+
I
R
.


## 3.3. Kinematic and Dynamic Models of Airship

The airship’s motion is subject to the net buoyancyB of airship and the atmosphere wind v
w. It can be assumed that the horizontal velocity of airship is equal to the atmospheric wind [17]. Therefore, the horizontal motion of airship can be assumed at the wind speed. An analysis of the vertical force balance on the airship is shown in Figure 2.Figure 2
Analysis of the vertical force balance on the airship.The thrust of airshipT is generated by vector thrust mechanism.The net buoyancy of airship is(20)
B
=
F
b
-
G
.The drag force of airshipD can be calculated by(21)
D
=
1
2
ρ
a
t
m
h
v
v
S
C
d
,where S
=
V
2
/
3 is the reference area and C
d is the drag coefficient.For the range of Reynolds numbers Re up to 107, the relationship between the drag coefficient C
d and Reynolds numbers Re was reported by Almedeij [18]. When Reynolds numbers are greater than 106, the drag coefficient is stable at the value of about 0.1, while for Reynolds numbers lower than 105 the constant drag coefficient of 0.5 can be used. A drag coefficient model for the ascent balloon in a Reynolds number range of 105~106 is explored by Conner and Arena [19], which can be expressed as follows:(22)
C
d
=
0.72
-
2.57
×
1
0
-
6
R
e
+
4.71
×
1
0
-
12
R
e
2
-
4.04
×
1
0
-
18
R
e
3
+
1.31
×
1
0
-
24
R
e
.The basic equations of airship during ascending can be expressed as follows:(23)
h
˙
=
v
,
v
˙
=
B
-
D
+
T
M
,
M
˙
=
m
˙
a
i
r
,where m
˙
a
i
r is the change rate of internal air mass. It is assumed that the airship only exhausts air and does not inhale air in this paper, so m
˙
a
i
r
≤
0.


## 3.4. Thermal Effects of Airship

The thermal effects of airship come from external and internal environment. External environment includes direct solar radiation, diffuse solar radiation, reflected radiation, infrared radiation, and forced convection between atmosphere and film. Internal environment includes infrared radiation among inner surface area elements and natural convection between inner gas and film. Thermal environment of airship is shown in Figure3.Figure 3
Thermal environment of airship.


### 3.4.1. Membrane Temperature Model

The surface of airship is divided intoM
×
N grids. The transient energy equation for each grid A
i
j can be expressed as follows [20]:(24)
d
s
ρ
s
c
T
˙
i
j
=
k
a
t
m
,
i
j
T
a
t
m
-
T
i
j
+
k
i
n
,
i
j
T
i
n
-
T
i
j
-
ϵ
1
+
ϵ
2
σ
T
i
j
4
+
∑
k
=
1
4
q
k
,
i
j
+
∑
I
=
1
M
∑
J
=
1
N
ϵ
2
E
I
J
A
i
j
·
∫
A
i
j
∫
A
I
J
cos
⁡
α
i
j
cos
⁡
α
I
J
π
l
2
d
A
i
j
d
A
I
J
,where d
s is the thickness of the film, ρ and c are the density and specific heat of the film, respectively, T
i
j is the temperature of element A
i
j, T
i
n is the inner gas temperature which contacts to the film, k
a
t
m
,
i
j is the heat transfer coefficient between atmosphere and film, k
i
n
,
i
j is the heat transfer coefficient between inner gas and film, ϵ
1 and ϵ
2 are the emissivity of the internal and external of film, respectively, E
I
J is the radiation heat flux of element A
I
J, α
i
j and α
I
J are the included angle between the internal normal vector n
→
i
j, n
→
I
J, and the connection of two elements, respectively, and q
k
,
i
j is the absorbed radiation heat flux of external film, including direct solar radiation q
1
,
i
j, diffuse solar radiation q
2
,
i
j, reflected radiation q
3
,
i
j, and infrared radiation q
4
,
i
j. The absorbed solar radiation [21] and absorbed infrared radiation [13] can be expressed as follows:(i)
The absorbed direct solar radiationq
1
,
i
j can be expressed as(25)
q
1
,
i
j
=
α
s
A
i
j
I
D
cos
⁡
β
i
j
,
whereα
s is the solar absorptivity of the film and β
i
j is the included angle between the element normal and the solar irradiation.
(ii)
The absorbed diffuse solar radiationq
2
,
i
j can be expressed as(26)
q
2
,
i
j
=
α
s
A
i
j
I
S
0.5
-
0.5
cos
⁡
ϕ
i
j
,
whereϕ
i
j is the included angle between the element normal and the gravity direction.
(iii)
The absorbed reflected solar radiationq
3
,
i
j can be expressed as(27)
q
3
,
i
j
=
α
s
A
i
j
I
R
0.5
+
0.5
cos
⁡
ϕ
i
j
.
(iv)
The absorbed infrared radiationq
4
,
i
j can be expressed as(28)
q
4
,
i
j
=
ϵ
2
A
i
j
σ
T
a
t
m
4
-
T
i
j
4
,
whereσ is the Stefan-Boltzmann constant.


### 3.4.2. Energy Equations of Inner Gas

Supposing that the temperature in one balloon is the same, according to the First Law of Thermodynamics, there are(29)
C
p
,
a
i
r
m
a
i
r
T
˙
a
i
r
=
Q
a
i
r
+
V
a
i
r
P
˙
i
n
-
C
p
,
a
i
r
T
a
i
r
m
˙
a
i
r
,
C
p
,
H
e
m
H
e
T
˙
H
e
=
Q
H
e
+
V
H
e
P
˙
i
n
,where C
p
,
a
i
r and C
p
,
H
e are the isobaric specific heats of air and helium, respectively. T
a
i
r and T
H
e are the temperature of air and helium, respectively. Q
a
i
r and Q
H
e are the quantity of heat convection of air and helium, respectively, which are given by(30)
Q
a
i
r
=
∑
i
=
1
M
∑
j
=
1
N
k
a
i
r
,
i
j
T
i
j
-
T
a
i
r
A
i
j
+
k
a
i
r
,
d
i
a
T
d
i
a
-
T
a
i
r
S
d
i
a
,
a
i
r
,
Q
H
e
=
∑
i
=
1
M
∑
j
=
1
N
k
H
e
,
i
j
T
i
j
-
T
H
e
A
i
j
+
k
H
e
,
d
i
a
T
d
i
a
-
T
H
e
S
d
i
a
,
H
e
,where T
d
i
a is the temperature of diaphragm between helium bag and air ballonet. k
a
i
r
,
i
j and k
a
i
r
,
d
i
a are the convection heat transfer coefficient between and film, diaphragm, and internal air. k
a
i
r
,
i
j, k
a
i
r
,
d
i
a, k
H
e
,
i
j and k
H
e
,
d
i
a are the convection heat transfer coefficients between internal air and film, internal air and diaphragm, helium and film, helium and diaphragm respectively. S
a
i
r
,
d
i
a and S
H
e
,
d
i
a are the contact areas between internal air and diaphragm, helium and diaphragm respectively. P
˙
i
n is the change rate of internal pressure which is given by(31)
V
P
˙
i
n
=
m
a
i
r
R
a
i
r
T
˙
a
i
r
+
m
H
e
R
H
e
T
˙
H
e
-
m
˙
a
i
r
R
a
i
r
T
a
i
r
.Combining (29)~(31), the change rate of inner air, helium, and internal pressure can be expressed as(32)
T
˙
a
i
r
=
Q
a
i
r
R
H
e
V
H
e
-
Q
H
e
R
H
e
V
a
i
r
-
C
p
,
H
e
Q
a
i
r
V
+
C
p
,
H
e
R
a
i
r
T
a
i
r
V
a
i
r
m
˙
a
i
r
-
C
p
,
a
i
r
R
H
e
T
a
i
r
V
H
e
m
˙
a
i
r
+
C
p
,
H
e
C
p
,
a
i
r
T
a
i
r
V
m
˙
a
i
r
m
a
i
r
C
p
,
H
e
R
a
i
r
V
a
i
r
+
C
p
,
a
i
r
R
H
e
V
H
e
-
C
p
,
H
e
C
p
,
a
i
r
V
,
T
˙
H
e
=
Q
H
e
R
a
i
r
V
a
i
r
-
Q
a
i
r
R
a
i
r
V
H
e
-
C
p
,
a
i
r
Q
H
e
V
+
2
C
p
,
a
i
r
R
a
i
r
T
a
i
r
V
H
e
m
˙
a
i
r
m
H
e
C
p
,
H
e
R
a
i
r
V
a
i
r
+
C
p
,
a
i
r
R
H
e
V
H
e
-
C
p
,
H
e
C
p
,
a
i
r
V
,
P
˙
i
n
=
-
C
p
,
H
e
Q
a
i
r
R
a
i
r
+
C
p
,
a
i
r
Q
H
e
R
H
e
-
2
C
p
,
H
e
C
p
,
a
i
r
R
a
i
r
T
a
i
r
m
˙
a
i
r
C
p
,
H
e
R
a
i
r
V
a
i
r
+
C
p
,
a
i
r
R
H
e
V
H
e
-
C
p
,
H
e
C
p
,
a
i
r
V
.


## 3.4.1. Membrane Temperature Model

The surface of airship is divided intoM
×
N grids. The transient energy equation for each grid A
i
j can be expressed as follows [20]:(24)
d
s
ρ
s
c
T
˙
i
j
=
k
a
t
m
,
i
j
T
a
t
m
-
T
i
j
+
k
i
n
,
i
j
T
i
n
-
T
i
j
-
ϵ
1
+
ϵ
2
σ
T
i
j
4
+
∑
k
=
1
4
q
k
,
i
j
+
∑
I
=
1
M
∑
J
=
1
N
ϵ
2
E
I
J
A
i
j
·
∫
A
i
j
∫
A
I
J
cos
⁡
α
i
j
cos
⁡
α
I
J
π
l
2
d
A
i
j
d
A
I
J
,where d
s is the thickness of the film, ρ and c are the density and specific heat of the film, respectively, T
i
j is the temperature of element A
i
j, T
i
n is the inner gas temperature which contacts to the film, k
a
t
m
,
i
j is the heat transfer coefficient between atmosphere and film, k
i
n
,
i
j is the heat transfer coefficient between inner gas and film, ϵ
1 and ϵ
2 are the emissivity of the internal and external of film, respectively, E
I
J is the radiation heat flux of element A
I
J, α
i
j and α
I
J are the included angle between the internal normal vector n
→
i
j, n
→
I
J, and the connection of two elements, respectively, and q
k
,
i
j is the absorbed radiation heat flux of external film, including direct solar radiation q
1
,
i
j, diffuse solar radiation q
2
,
i
j, reflected radiation q
3
,
i
j, and infrared radiation q
4
,
i
j. The absorbed solar radiation [21] and absorbed infrared radiation [13] can be expressed as follows:(i)
The absorbed direct solar radiationq
1
,
i
j can be expressed as(25)
q
1
,
i
j
=
α
s
A
i
j
I
D
cos
⁡
β
i
j
,
whereα
s is the solar absorptivity of the film and β
i
j is the included angle between the element normal and the solar irradiation.
(ii)
The absorbed diffuse solar radiationq
2
,
i
j can be expressed as(26)
q
2
,
i
j
=
α
s
A
i
j
I
S
0.5
-
0.5
cos
⁡
ϕ
i
j
,
whereϕ
i
j is the included angle between the element normal and the gravity direction.
(iii)
The absorbed reflected solar radiationq
3
,
i
j can be expressed as(27)
q
3
,
i
j
=
α
s
A
i
j
I
R
0.5
+
0.5
cos
⁡
ϕ
i
j
.
(iv)
The absorbed infrared radiationq
4
,
i
j can be expressed as(28)
q
4
,
i
j
=
ϵ
2
A
i
j
σ
T
a
t
m
4
-
T
i
j
4
,
whereσ is the Stefan-Boltzmann constant.


## 3.4.2. Energy Equations of Inner Gas

Supposing that the temperature in one balloon is the same, according to the First Law of Thermodynamics, there are(29)
C
p
,
a
i
r
m
a
i
r
T
˙
a
i
r
=
Q
a
i
r
+
V
a
i
r
P
˙
i
n
-
C
p
,
a
i
r
T
a
i
r
m
˙
a
i
r
,
C
p
,
H
e
m
H
e
T
˙
H
e
=
Q
H
e
+
V
H
e
P
˙
i
n
,where C
p
,
a
i
r and C
p
,
H
e are the isobaric specific heats of air and helium, respectively. T
a
i
r and T
H
e are the temperature of air and helium, respectively. Q
a
i
r and Q
H
e are the quantity of heat convection of air and helium, respectively, which are given by(30)
Q
a
i
r
=
∑
i
=
1
M
∑
j
=
1
N
k
a
i
r
,
i
j
T
i
j
-
T
a
i
r
A
i
j
+
k
a
i
r
,
d
i
a
T
d
i
a
-
T
a
i
r
S
d
i
a
,
a
i
r
,
Q
H
e
=
∑
i
=
1
M
∑
j
=
1
N
k
H
e
,
i
j
T
i
j
-
T
H
e
A
i
j
+
k
H
e
,
d
i
a
T
d
i
a
-
T
H
e
S
d
i
a
,
H
e
,where T
d
i
a is the temperature of diaphragm between helium bag and air ballonet. k
a
i
r
,
i
j and k
a
i
r
,
d
i
a are the convection heat transfer coefficient between and film, diaphragm, and internal air. k
a
i
r
,
i
j, k
a
i
r
,
d
i
a, k
H
e
,
i
j and k
H
e
,
d
i
a are the convection heat transfer coefficients between internal air and film, internal air and diaphragm, helium and film, helium and diaphragm respectively. S
a
i
r
,
d
i
a and S
H
e
,
d
i
a are the contact areas between internal air and diaphragm, helium and diaphragm respectively. P
˙
i
n is the change rate of internal pressure which is given by(31)
V
P
˙
i
n
=
m
a
i
r
R
a
i
r
T
˙
a
i
r
+
m
H
e
R
H
e
T
˙
H
e
-
m
˙
a
i
r
R
a
i
r
T
a
i
r
.Combining (29)~(31), the change rate of inner air, helium, and internal pressure can be expressed as(32)
T
˙
a
i
r
=
Q
a
i
r
R
H
e
V
H
e
-
Q
H
e
R
H
e
V
a
i
r
-
C
p
,
H
e
Q
a
i
r
V
+
C
p
,
H
e
R
a
i
r
T
a
i
r
V
a
i
r
m
˙
a
i
r
-
C
p
,
a
i
r
R
H
e
T
a
i
r
V
H
e
m
˙
a
i
r
+
C
p
,
H
e
C
p
,
a
i
r
T
a
i
r
V
m
˙
a
i
r
m
a
i
r
C
p
,
H
e
R
a
i
r
V
a
i
r
+
C
p
,
a
i
r
R
H
e
V
H
e
-
C
p
,
H
e
C
p
,
a
i
r
V
,
T
˙
H
e
=
Q
H
e
R
a
i
r
V
a
i
r
-
Q
a
i
r
R
a
i
r
V
H
e
-
C
p
,
a
i
r
Q
H
e
V
+
2
C
p
,
a
i
r
R
a
i
r
T
a
i
r
V
H
e
m
˙
a
i
r
m
H
e
C
p
,
H
e
R
a
i
r
V
a
i
r
+
C
p
,
a
i
r
R
H
e
V
H
e
-
C
p
,
H
e
C
p
,
a
i
r
V
,
P
˙
i
n
=
-
C
p
,
H
e
Q
a
i
r
R
a
i
r
+
C
p
,
a
i
r
Q
H
e
R
H
e
-
2
C
p
,
H
e
C
p
,
a
i
r
R
a
i
r
T
a
i
r
m
˙
a
i
r
C
p
,
H
e
R
a
i
r
V
a
i
r
+
C
p
,
a
i
r
R
H
e
V
H
e
-
C
p
,
H
e
C
p
,
a
i
r
V
.


## 3.5. Actuator Model

Actuator includes exhaust valve and vector thrust mechanism. The exhaust valve expels inner air to maintain the pressure difference while vector thrust mechanism generates the vertical thrust to help control rising velocity.The parameters of exhaust valve are as follows: the diameter of valve isd
a
i
r, the adjustable opening area is S
a
i
r
=
π
d
a
i
r, the opening degree is f
a
i
r
(
t
), and flow resistance is C
a
i
r. According to Bernoulli’s principle, there is(33)
m
˙
a
i
r
t
=
2
P
i
n
t
Δ
P
t
R
a
i
r
T
a
i
r
t
S
a
i
r
C
a
i
r
f
a
i
r
t
.The effect of exhausting air is relevant not only to opening degree but also to the current thermodynamic states of inner gas.The vector thrust mechanism can continuously generate the vertical thrustT, and the range of T is 0~1000 N.


## 4. Realization of Coordinated Control

Based on the above models, the coordinated control of pressure difference and rising velocity for stratospheric airship is studied by addressing the following two issues:(1)
consider the interference caused by solar radiation and other thermal processes;
(2)
the coordinated control of pressure difference and rising velocity can be achieved by tracking the corresponding reference trajectories smoothly, thereby ensuring that the airship rises from sea level to the altitude of 20 km steadily and rapidly.


### 4.1. Design of Reference Trajectories

#### 4.1.1. Pressure Difference Reference Trajectory

Assuming that the initial pressure difference is 500 Pa and the expected pressure difference is 350 Pa, so the reference trajectory of pressure difference can be designed as follows:(34)
r
Δ
P
t
=
350
+
50
e
-
0.05
t
+
100
e
-
0.01
t
.


#### 4.1.2. Rising Velocity Reference Trajectory

The wind speed will increase obviously with height above the 7 km, especially in the strong wind zone 11~12 km. As a low-velocity aircraft, the stratospheric airship has a high sensibility to wind disturbance, so it should pass through this harsh environment quickly in the limited time. The reference trajectory of rising velocity is designed as follows.(i)
The reference trajectory of acceleration and uniform velocity is as follows:
when0
≤
h
≤
7000
m,(35)
r
v
t
=
v
r
1
+
e
-
0.02
t
-
2
e
-
0.01
t
,
and when7000
m
<
h
≤
20000
m,(36)
r
v
t
=
v
r
2
+
e
-
0.01
t
-
t
s
1
-
2
e
-
0.005
t
-
t
s
1
.
(ii)
The reference trajectory of deceleration velocity is(37)
r
v
t
=
2
v
r
2
e
-
0.01
t
-
t
s
2
-
e
-
0.02
t
-
t
s
2
,
wherev
r
=
0.5
m
/
s. t
s
1 and t
s
2 are switching times.


#### 4.1.3. Switching Time

The switching time ist
s
1 when airship reaches 7 km. Because the rising displacement h
d is constant during the process of decelerating, the switching time is t
s
2 when airship reaches (
20
k
m
-
h
d
).


### 4.2. Coordinated Control of Airship

Based on the the current internal and external pressure, environment control system (ECS) controls exhaust or inspiratory valve to regulate the pressure difference. The airship geometry is substantially maintained by remaining a certain pressure difference. For tracking the reference trajectory of rising velocity, flight control system (FCS) controls vector thrust mechanism to generate the vertical thrust for helping control rising velocity. The logic diagram of coordinated control of pressure difference and rising velocity is shown in Figure4. Because exhausting or inhaling air will change the airship mass to affect rising velocity and the rising velocity will affect pressure difference through changing height, these two processes are highly coupled. The controllers of pressure difference and rising velocity are designed for tracking the reference trajectories by the incomplete differential control algorithm and pole placement method, respectively, so the rules-based changes of pressure difference and rising velocity are conducive to decouple and realize the coordinated control of pressure difference and rising velocity.Figure 4
Logic diagram of coordinated control.


#### 4.2.1. Control of Pressure Difference

The error between the actual pressure difference and the pressure difference reference trajectory is(38)
e
Δ
P
t
=
Δ
P
t
-
r
Δ
P
t
.Based on analog PID control algorithm, the continuous time t is approximated by a series of sampling points k
T, T is the sampling period, k is the sampling sequence number, k
=
1,2
,
…, and the differential is approximated by first order backward difference, so there is(39)
t
≈
k
T
;
d
e
Δ
P
t
d
t
≈
e
Δ
P
k
T
-
e
Δ
P
k
-
1
T
T
=
e
Δ
P
k
-
e
Δ
P
k
-
1
T
.According to the feedback quantity of e
Δ
P
(
t
) and the change rate of e
Δ
P
(
t
), the PD controller is designed in discrete domain as follows:(40)
u
′
k
=
k
p
e
Δ
P
k
+
k
d
e
Δ
P
k
-
e
Δ
P
k
-
1
T
,where k
p is the proportional coefficient and k
d is the differential coefficient. e
Δ
P
(
k
) and e
Δ
P
(
k
-
1
) are error signal of pressure difference at k sampling point and at (
k
-
1
) sampling point, respectively. The derivative action of PD controller can improve the dynamic characteristics of the system, but it is also particularly sensitive to high-frequency signal interference. Therefore, the incomplete differential PD controller is designed by using a low-pass filter to achieve the desired filtering effect. The transfer function of low-pass filter is(41)
D
f
s
=
1
1
+
T
f
s
,where T
f
=
k
d
/
k
p
/
N is the filtering time constant and the range of N is 2~20.According to Figure5, there is(42)
T
f
d
u
t
d
t
+
u
t
=
u
′
t
.The incomplete differential PD controller is designed by discretizing (42):(43)
u
k
=
T
T
f
+
T
u
k
-
1
+
T
f
T
f
+
T
u
′
k
.So the change rate of m
a
i
r is(44)
m
˙
a
i
r
t
=
u
t
,
e
Δ
P
t
>
0
;
0
,
e
Δ
P
t
≤
0
.Figure 5
Structure diagram of incomplete differential PD controller.According to (33), the opening degree of exhaust valve is(45)
f
a
i
r
t
=
m
˙
a
i
r
t
2
P
i
n
t
Δ
P
t
/
R
a
i
r
T
a
i
r
t
S
a
i
r
C
a
i
r
.


#### 4.2.2. Control of Rising Velocity

The error between current rising velocity and the rising velocity reference trajectory is(46)
e
v
t
=
v
t
-
r
v
t
.So the change rate of e
v can be expressed as(47)
e
˙
v
t
=
v
˙
t
-
r
˙
v
t
.Combining (23) and (47), there is(48)
e
˙
v
t
=
B
t
-
D
t
+
T
t
M
t
-
r
˙
v
t
.Based on pole placement method, e
˙
v can also be expressed as(49)
e
˙
v
t
=
-
K
1
e
v
t
,where K
1 is a constant and K
1
>
0. And(50)
e
v
t
=
e
v
0
e
-
K
1
t
,
l
i
m
t
→
∞
e
v
t
=
0
.Combining (48) and (49), the vertical thrust is(51)
T
t
=
M
t
-
K
1
e
v
t
+
r
˙
v
t
-
B
t
-
D
t
,
e
v
t
<
0
;
0
,
e
v
t
≥
0
.


## 4.1. Design of Reference Trajectories

### 4.1.1. Pressure Difference Reference Trajectory

Assuming that the initial pressure difference is 500 Pa and the expected pressure difference is 350 Pa, so the reference trajectory of pressure difference can be designed as follows:(34)
r
Δ
P
t
=
350
+
50
e
-
0.05
t
+
100
e
-
0.01
t
.


### 4.1.2. Rising Velocity Reference Trajectory

The wind speed will increase obviously with height above the 7 km, especially in the strong wind zone 11~12 km. As a low-velocity aircraft, the stratospheric airship has a high sensibility to wind disturbance, so it should pass through this harsh environment quickly in the limited time. The reference trajectory of rising velocity is designed as follows.(i)
The reference trajectory of acceleration and uniform velocity is as follows:
when0
≤
h
≤
7000
m,(35)
r
v
t
=
v
r
1
+
e
-
0.02
t
-
2
e
-
0.01
t
,
and when7000
m
<
h
≤
20000
m,(36)
r
v
t
=
v
r
2
+
e
-
0.01
t
-
t
s
1
-
2
e
-
0.005
t
-
t
s
1
.
(ii)
The reference trajectory of deceleration velocity is(37)
r
v
t
=
2
v
r
2
e
-
0.01
t
-
t
s
2
-
e
-
0.02
t
-
t
s
2
,
wherev
r
=
0.5
m
/
s. t
s
1 and t
s
2 are switching times.


### 4.1.3. Switching Time

The switching time ist
s
1 when airship reaches 7 km. Because the rising displacement h
d is constant during the process of decelerating, the switching time is t
s
2 when airship reaches (
20
k
m
-
h
d
).


## 4.1.1. Pressure Difference Reference Trajectory

Assuming that the initial pressure difference is 500 Pa and the expected pressure difference is 350 Pa, so the reference trajectory of pressure difference can be designed as follows:(34)
r
Δ
P
t
=
350
+
50
e
-
0.05
t
+
100
e
-
0.01
t
.


## 4.1.2. Rising Velocity Reference Trajectory

The wind speed will increase obviously with height above the 7 km, especially in the strong wind zone 11~12 km. As a low-velocity aircraft, the stratospheric airship has a high sensibility to wind disturbance, so it should pass through this harsh environment quickly in the limited time. The reference trajectory of rising velocity is designed as follows.(i)
The reference trajectory of acceleration and uniform velocity is as follows:
when0
≤
h
≤
7000
m,(35)
r
v
t
=
v
r
1
+
e
-
0.02
t
-
2
e
-
0.01
t
,
and when7000
m
<
h
≤
20000
m,(36)
r
v
t
=
v
r
2
+
e
-
0.01
t
-
t
s
1
-
2
e
-
0.005
t
-
t
s
1
.
(ii)
The reference trajectory of deceleration velocity is(37)
r
v
t
=
2
v
r
2
e
-
0.01
t
-
t
s
2
-
e
-
0.02
t
-
t
s
2
,
wherev
r
=
0.5
m
/
s. t
s
1 and t
s
2 are switching times.


## 4.1.3. Switching Time

The switching time ist
s
1 when airship reaches 7 km. Because the rising displacement h
d is constant during the process of decelerating, the switching time is t
s
2 when airship reaches (
20
k
m
-
h
d
).


## 4.2. Coordinated Control of Airship

Based on the the current internal and external pressure, environment control system (ECS) controls exhaust or inspiratory valve to regulate the pressure difference. The airship geometry is substantially maintained by remaining a certain pressure difference. For tracking the reference trajectory of rising velocity, flight control system (FCS) controls vector thrust mechanism to generate the vertical thrust for helping control rising velocity. The logic diagram of coordinated control of pressure difference and rising velocity is shown in Figure4. Because exhausting or inhaling air will change the airship mass to affect rising velocity and the rising velocity will affect pressure difference through changing height, these two processes are highly coupled. The controllers of pressure difference and rising velocity are designed for tracking the reference trajectories by the incomplete differential control algorithm and pole placement method, respectively, so the rules-based changes of pressure difference and rising velocity are conducive to decouple and realize the coordinated control of pressure difference and rising velocity.Figure 4
Logic diagram of coordinated control.


### 4.2.1. Control of Pressure Difference

The error between the actual pressure difference and the pressure difference reference trajectory is(38)
e
Δ
P
t
=
Δ
P
t
-
r
Δ
P
t
.Based on analog PID control algorithm, the continuous time t is approximated by a series of sampling points k
T, T is the sampling period, k is the sampling sequence number, k
=
1,2
,
…, and the differential is approximated by first order backward difference, so there is(39)
t
≈
k
T
;
d
e
Δ
P
t
d
t
≈
e
Δ
P
k
T
-
e
Δ
P
k
-
1
T
T
=
e
Δ
P
k
-
e
Δ
P
k
-
1
T
.According to the feedback quantity of e
Δ
P
(
t
) and the change rate of e
Δ
P
(
t
), the PD controller is designed in discrete domain as follows:(40)
u
′
k
=
k
p
e
Δ
P
k
+
k
d
e
Δ
P
k
-
e
Δ
P
k
-
1
T
,where k
p is the proportional coefficient and k
d is the differential coefficient. e
Δ
P
(
k
) and e
Δ
P
(
k
-
1
) are error signal of pressure difference at k sampling point and at (
k
-
1
) sampling point, respectively. The derivative action of PD controller can improve the dynamic characteristics of the system, but it is also particularly sensitive to high-frequency signal interference. Therefore, the incomplete differential PD controller is designed by using a low-pass filter to achieve the desired filtering effect. The transfer function of low-pass filter is(41)
D
f
s
=
1
1
+
T
f
s
,where T
f
=
k
d
/
k
p
/
N is the filtering time constant and the range of N is 2~20.According to Figure5, there is(42)
T
f
d
u
t
d
t
+
u
t
=
u
′
t
.The incomplete differential PD controller is designed by discretizing (42):(43)
u
k
=
T
T
f
+
T
u
k
-
1
+
T
f
T
f
+
T
u
′
k
.So the change rate of m
a
i
r is(44)
m
˙
a
i
r
t
=
u
t
,
e
Δ
P
t
>
0
;
0
,
e
Δ
P
t
≤
0
.Figure 5
Structure diagram of incomplete differential PD controller.According to (33), the opening degree of exhaust valve is(45)
f
a
i
r
t
=
m
˙
a
i
r
t
2
P
i
n
t
Δ
P
t
/
R
a
i
r
T
a
i
r
t
S
a
i
r
C
a
i
r
.


### 4.2.2. Control of Rising Velocity

The error between current rising velocity and the rising velocity reference trajectory is(46)
e
v
t
=
v
t
-
r
v
t
.So the change rate of e
v can be expressed as(47)
e
˙
v
t
=
v
˙
t
-
r
˙
v
t
.Combining (23) and (47), there is(48)
e
˙
v
t
=
B
t
-
D
t
+
T
t
M
t
-
r
˙
v
t
.Based on pole placement method, e
˙
v can also be expressed as(49)
e
˙
v
t
=
-
K
1
e
v
t
,where K
1 is a constant and K
1
>
0. And(50)
e
v
t
=
e
v
0
e
-
K
1
t
,
l
i
m
t
→
∞
e
v
t
=
0
.Combining (48) and (49), the vertical thrust is(51)
T
t
=
M
t
-
K
1
e
v
t
+
r
˙
v
t
-
B
t
-
D
t
,
e
v
t
<
0
;
0
,
e
v
t
≥
0
.


## 4.2.1. Control of Pressure Difference

The error between the actual pressure difference and the pressure difference reference trajectory is(38)
e
Δ
P
t
=
Δ
P
t
-
r
Δ
P
t
.Based on analog PID control algorithm, the continuous time t is approximated by a series of sampling points k
T, T is the sampling period, k is the sampling sequence number, k
=
1,2
,
…, and the differential is approximated by first order backward difference, so there is(39)
t
≈
k
T
;
d
e
Δ
P
t
d
t
≈
e
Δ
P
k
T
-
e
Δ
P
k
-
1
T
T
=
e
Δ
P
k
-
e
Δ
P
k
-
1
T
.According to the feedback quantity of e
Δ
P
(
t
) and the change rate of e
Δ
P
(
t
), the PD controller is designed in discrete domain as follows:(40)
u
′
k
=
k
p
e
Δ
P
k
+
k
d
e
Δ
P
k
-
e
Δ
P
k
-
1
T
,where k
p is the proportional coefficient and k
d is the differential coefficient. e
Δ
P
(
k
) and e
Δ
P
(
k
-
1
) are error signal of pressure difference at k sampling point and at (
k
-
1
) sampling point, respectively. The derivative action of PD controller can improve the dynamic characteristics of the system, but it is also particularly sensitive to high-frequency signal interference. Therefore, the incomplete differential PD controller is designed by using a low-pass filter to achieve the desired filtering effect. The transfer function of low-pass filter is(41)
D
f
s
=
1
1
+
T
f
s
,where T
f
=
k
d
/
k
p
/
N is the filtering time constant and the range of N is 2~20.According to Figure5, there is(42)
T
f
d
u
t
d
t
+
u
t
=
u
′
t
.The incomplete differential PD controller is designed by discretizing (42):(43)
u
k
=
T
T
f
+
T
u
k
-
1
+
T
f
T
f
+
T
u
′
k
.So the change rate of m
a
i
r is(44)
m
˙
a
i
r
t
=
u
t
,
e
Δ
P
t
>
0
;
0
,
e
Δ
P
t
≤
0
.Figure 5
Structure diagram of incomplete differential PD controller.According to (33), the opening degree of exhaust valve is(45)
f
a
i
r
t
=
m
˙
a
i
r
t
2
P
i
n
t
Δ
P
t
/
R
a
i
r
T
a
i
r
t
S
a
i
r
C
a
i
r
.


## 4.2.2. Control of Rising Velocity

The error between current rising velocity and the rising velocity reference trajectory is(46)
e
v
t
=
v
t
-
r
v
t
.So the change rate of e
v can be expressed as(47)
e
˙
v
t
=
v
˙
t
-
r
˙
v
t
.Combining (23) and (47), there is(48)
e
˙
v
t
=
B
t
-
D
t
+
T
t
M
t
-
r
˙
v
t
.Based on pole placement method, e
˙
v can also be expressed as(49)
e
˙
v
t
=
-
K
1
e
v
t
,where K
1 is a constant and K
1
>
0. And(50)
e
v
t
=
e
v
0
e
-
K
1
t
,
l
i
m
t
→
∞
e
v
t
=
0
.Combining (48) and (49), the vertical thrust is(51)
T
t
=
M
t
-
K
1
e
v
t
+
r
˙
v
t
-
B
t
-
D
t
,
e
v
t
<
0
;
0
,
e
v
t
≥
0
.


## 5. Numerical Results and Discussions

### 5.1. Simulation Parameters

The classical examples are described by [22, 23], so simulation parameters are considered as follows. The take-off point of airship is Xiamen, China, being confined by longitudes 11
8
∘
04
′
04
′
′ and latitudes 2
4
∘
26
′
46
′
′, and the take-off time is 8:00 am on the autumnal equinox. A double ellipsoid model of airship with the same minor axis b
=
3.75
m, the long axis lengths are a
1
=
8.33
m and a
2
=
16.66
m, respectively, the volume is V
=
736
m
3, the initial total mass is m
(
0
)
=
900
k
g, the helium mass is m
H
e
=
9
k
g, the initial air mass is m
a
i
r
(
0
)
=
839
k
g, the structure mass is m
s
=
52
k
g, the gas constants of air and helium are R
a
i
r
=
287 and R
H
e
=
2078, the isobaric specific heats of air and helium are C
p
,
a
i
r
=
1011
J
/
(
k
g
·
K
) and C
p
,
H
e
=
5225
J
/
(
k
g
·
K
), the density of film ρ
s
k
i
n
=
0.44
k
g
/
m
3, the solar absorption coefficient of film α
s is 0.13, and the emissivity of the internal and external of flim is ϵ
1
=
ϵ
2
=
0.43.


### 5.2. Results and Discussions

Figure6 shows that airship rises to the altitude of 20 km smoothly within 8 hours and then keeps itself at this height. The profiles of rising velocity and pressure difference are presented in Figures 7~9. The velocity tracks the reference trajectory well, while the pressure difference is maintained in the safe range. The mass change and the opening degree of exhaust valve are presented in Figures 10~12. Figures 13 and 14 show the profiles of vector thrust and power of actuator, respectively. Figures 15 and 16 illustrate the profile of temperature difference between the inner gas and atmosphere. Because the thermal models, such as solar radiation and membrane temperature change, are not used to design the controllers, these results show that the control scheme not only can achieve the coordinated control of pressure difference and rising velocity but also has good robustness against the thermodynamic disturbance.Figure 6
Flight height profile during airship ascending.Figure 7
Rising velocity profile during airship ascending.Figure 8
Internal and external pressure profile during airship ascending.Figure 9
Pressure difference profile during airship ascending.Figure 10
Inner air mass profile during airship ascending.It is noticeable that Figures15 and 16 show that the temperature of the inner air is obviously less than the temperature of helium, while the temperature of helium is less than the temperature of the surrounding atmosphere. This phenomenon is called “supercool” in airship research. The reason is that to track the velocity reference trajectory and guarantee pressure difference during ascending the airship must exhaust inner air continuously. With thermodynamic theory, exhausting air leads to the decrease of inner air’s temperature and pressure. Meanwhile, helium also expands and cools. Comparing Figures 11 and 16, we can find that the profile of the exhausting rate of inner air is highly consistent with the profile of the temperature difference.Figure 11
The change rate of inner air mass profile during airship ascending.Figure 12
The opening degree of exhaust valve profile during airship ascending.Figure 13
Vector thrust of actuator profile during airship ascending.Figure 14
Power of actuator profile during airship ascending.Figure 15
Temperature profile of inner air, helium, and atmosphere during airship ascending.Figure 16
The profile of temperature difference between inner air, helium, and atmosphere during airship ascending.“Supercool” is highly relative to the rising velocity. It is also affected by factors such as solar radiation and convection heat transfer. Generally, the faster the airship ascends, the more quickly the inner air’s temperature drops. So, it is important to design the rising velocity reference trajectory to overcome “supercool” and guarantee the temperature difference between the inside and outside of airship within a reasonable range.Opposite to “supercool,” “superheat” phenomenon is caused by inhaling air continuously and plays an important role in the airship’s descending.


## 5.1. Simulation Parameters

The classical examples are described by [22, 23], so simulation parameters are considered as follows. The take-off point of airship is Xiamen, China, being confined by longitudes 11
8
∘
04
′
04
′
′ and latitudes 2
4
∘
26
′
46
′
′, and the take-off time is 8:00 am on the autumnal equinox. A double ellipsoid model of airship with the same minor axis b
=
3.75
m, the long axis lengths are a
1
=
8.33
m and a
2
=
16.66
m, respectively, the volume is V
=
736
m
3, the initial total mass is m
(
0
)
=
900
k
g, the helium mass is m
H
e
=
9
k
g, the initial air mass is m
a
i
r
(
0
)
=
839
k
g, the structure mass is m
s
=
52
k
g, the gas constants of air and helium are R
a
i
r
=
287 and R
H
e
=
2078, the isobaric specific heats of air and helium are C
p
,
a
i
r
=
1011
J
/
(
k
g
·
K
) and C
p
,
H
e
=
5225
J
/
(
k
g
·
K
), the density of film ρ
s
k
i
n
=
0.44
k
g
/
m
3, the solar absorption coefficient of film α
s is 0.13, and the emissivity of the internal and external of flim is ϵ
1
=
ϵ
2
=
0.43.


## 5.2. Results and Discussions

Figure6 shows that airship rises to the altitude of 20 km smoothly within 8 hours and then keeps itself at this height. The profiles of rising velocity and pressure difference are presented in Figures 7~9. The velocity tracks the reference trajectory well, while the pressure difference is maintained in the safe range. The mass change and the opening degree of exhaust valve are presented in Figures 10~12. Figures 13 and 14 show the profiles of vector thrust and power of actuator, respectively. Figures 15 and 16 illustrate the profile of temperature difference between the inner gas and atmosphere. Because the thermal models, such as solar radiation and membrane temperature change, are not used to design the controllers, these results show that the control scheme not only can achieve the coordinated control of pressure difference and rising velocity but also has good robustness against the thermodynamic disturbance.Figure 6
Flight height profile during airship ascending.Figure 7
Rising velocity profile during airship ascending.Figure 8
Internal and external pressure profile during airship ascending.Figure 9
Pressure difference profile during airship ascending.Figure 10
Inner air mass profile during airship ascending.It is noticeable that Figures15 and 16 show that the temperature of the inner air is obviously less than the temperature of helium, while the temperature of helium is less than the temperature of the surrounding atmosphere. This phenomenon is called “supercool” in airship research. The reason is that to track the velocity reference trajectory and guarantee pressure difference during ascending the airship must exhaust inner air continuously. With thermodynamic theory, exhausting air leads to the decrease of inner air’s temperature and pressure. Meanwhile, helium also expands and cools. Comparing Figures 11 and 16, we can find that the profile of the exhausting rate of inner air is highly consistent with the profile of the temperature difference.Figure 11
The change rate of inner air mass profile during airship ascending.Figure 12
The opening degree of exhaust valve profile during airship ascending.Figure 13
Vector thrust of actuator profile during airship ascending.Figure 14
Power of actuator profile during airship ascending.Figure 15
Temperature profile of inner air, helium, and atmosphere during airship ascending.Figure 16
The profile of temperature difference between inner air, helium, and atmosphere during airship ascending.“Supercool” is highly relative to the rising velocity. It is also affected by factors such as solar radiation and convection heat transfer. Generally, the faster the airship ascends, the more quickly the inner air’s temperature drops. So, it is important to design the rising velocity reference trajectory to overcome “supercool” and guarantee the temperature difference between the inside and outside of airship within a reasonable range.Opposite to “supercool,” “superheat” phenomenon is caused by inhaling air continuously and plays an important role in the airship’s descending.


## 6. Conclusions

During stratospheric airship ascending, pressure difference and rising velocity are highly coupled. In this paper, a coordinated control strategy, in which the priority is given to control the pressure difference, is presented for the airship equipped with a vector thrust. The airship exhausts the internal air to maintain the pressure difference definitely in the safe range while the vector thrust generates the vertical thrust coworks with the static buoyancy to guarantee the rising velocity tracking a given reference velocity. In the control strategy, the thermal effects such as solar radiation are regarded as disturbances while mathematical models are established to evaluate these thermal effects. Simulations show that the coordinated control strategy can meet the requirements of pressure difference and vertical velocity during ascending and be robust enough against the thermal disturbance.



---

*Source: 923156-2015-06-11.xml*