643
1 INTRODUCTION
The article presents the mathematical analysis of the
pneumatic cascade and the pneumatic membrane
actuator described with the differential calculus of
non-integer orders (Fractional calculus) [1], [2], [3], [5],
[6], [7], [8], [9], [10], [13] and [15].
Differential equations of integer and non-integer
order were introduced and became the basis for
deriving equations describing time characteristics
(pulse and step) and frequency characteristics
(logarithmic amplitude and phase characteristics) for
each tested pneumatic system [11]. Next, simulations
of derived equations were performed using Microsoft
Excel and MATLAB & Simulink software, obtaining
time and frequency characteristics of the tested
systems for integer and non-integer orders [12], [14],
[16], [17] and [18].
2 MATHEMATICAL MODEL OF A PNEUMATIC
CASCADE DESCRIBED WITH A FRACTIONAL
CALCULUS
Fig. 1 shows the diagram of the analyzed pneumatic
cascade:
Figure 1. Diagram of a two-chamber pneumatic cascade [11]
Fig. 2 shows a block diagram of a two-chamber
pneumatic cascade:
Mathematical Models of a Pneumatic Cascade and
P
neumatic Membrane Actuator Described with a
F
ractional Calculus
M
. Luft, Z. Łukasik & D. Pietruszczak
Kazimierz Pulaski University of Technology and Humanities in Radom,
Radom, Poland
ABSTRACT: The paper presents the analysis of dynamic properties of pneumatic systems such as: cascade
connection of membrane pressure transmitters and a pneumatic membrane actuator by means of differential
equations of integer and non-integer order. The analyzed systems were described from the time perspective by
means of step response, and in terms of frequency with the help of the Bode plot, i.e. logarithmic magnitude
and phase responses.
Each response was determined using differential equations of non-integer order.
To determine the responses, the interactive Simulink package was an irreplaceable programming tool built on
the basis of the MATLAB program, which enables the analysis and synthesis of continuous dynamic systems.
http://www.transnav.eu
the
International Journal
on Marin
e Navigation
and Safety of Sea Transportation
Volume 14
Number 3
September 2020
DOI:
10.12716/1001.14.03.16
644
Figure 2. Block diagram of a two-chamber pneumatic
cascade [11]
Assuming the linearity of the model, the equation
describing the dynamics of the diaphragm pressure
transmitter can be written in the form of a system of
differential equations:
( )
( )
(
) (
)
2
11
22
11 1 1 1 0
2
2
dpt dpt
pt pt
dt
dt
ξω ω ω
+ +=
(1a)
( ) ( )
( ) (
)
2
22
22
22 2 21
2
2
dpt dpt
pt p t
dt
dt
ξω ω ω
+ +=
(1b)
where:
( )
1,2
rad
s
ω
- the pulsatance of another elementary
pneumatic system,
1,2
1
ξ
<
- the damping ratio of another elementary
pneumatic system included in the pneumatic cascade.
22
1,2
1,2
1,2 1,2
1, 2 1, 2
3
1
4
pp pp
rc
lV
LC
π
ω
= =
(2a)
1,2 1,2
2
1, 2 1, 2 1,2
1,2 1,2
1,2
2 22
1,2
1,2
3
3
22
2
pp pp
lV
LC
lV
rc
rc
η
ω
η
π
ζ
ρ
πρ
= = =
(2b)
wherein:
25
p
C Ns m
−
- pneumatic capacity in another element
of the pneumatic system;
31
p
L mN
−
- pneumatic induction in another element
of the pneumatic system;
5
p
R Nsm
−
- flow resistance in another element of
the pneumatic system;
3
Vm
- volume of another transducer chamber,
m
c
s
- speed of sound in the gas filling the system;
[ ]
lm
- length of another inlet pipe;
[ ]
rm
- radius of another inlet tube;
3
kgm
ρ
−
- gas density;
11
kgm s
η
−−
- dynamic viscosity.
Equations (1) written with the help of fractional
calculus take the following form:
( ) ( ) ( ) ( )
( ) ( ) ( ) ( )
2 22
0 1 110 1 1 1 1 0
2 22
0 2 2 20 2 2 2 2 1
2
2
RL v RL v
tt
RL v RL v
tt
Dpt Dpt pt pt
Dpt Dpt pt pt
ξω ω ω
ξω ω ω
+ +=
+ +=
(3)
where:
0v >
.
Using the Laplace transform to equation (3), for
zero initial conditions, we obtain:
( )
(
) ( )
(
)
( )
(
)
2 22
11 1 1 1 0
2 22
22 2 2 21
2
2
vv
vv
s s ps ps
s s ps ps
ξω ω ω
ξω ω ω
++ =
++ =
(4)
Thus the transfer function of non-integer order of
the analyzed pneumatic system is obtained:
( )
( )
(
)
(
)
(
)
( )
(
)
( )
2
1
1
1
22
0
11 1
2
2
2
2
22
1
22 2
2
2
v
vv
v
vv
ps
Gs
ps
ss
ps
Gs
ps
ss
ω
ξω ω
ω
ξω ω
= =
++
= =
++
(5)
The transfer function of the analyzed system takes
the form:
( )
( )
( )
( )
( )
( )
( )
( )
( )
( )
( )
( )
( )
2
12
0
22
12
4 3 2 22
1 1 2 2 1 12 1 2 2
2 2 22
112 212 12
22 4
22
v vv
v
vv v
v
ps
Gs GsGs
ps
Gs
ss s
s
ωω
ξω ξω ω ξξωω ω
ξωω ξωω ωω
= =
=
+ + ++ + +
++
(6)
For the formula (6), we obtain the spectral transfer
function of the tested transducer:
( )
( ) ( )
( )
( )
( ) ( )
( )
22
12
4
3
11 2 2
2 22
1 12 1 2 2
2 2 22
112 212 12
cos 2 sin 2
33
2 2 cos sin
22
4 cos sin
2 2 cos sin
22
v
v
v
v
v
Gj
vj v
vv
j
vj v
vv
j
ωω
ω
ωπ π
ππ
ξω ξω ω
ω ξξ ωω ω ω π π
ππ
ξωω ξωω ω ωω
=
++
+ ++
++ + +
++ + +
(7)
By making elementary transformations, the real
and imaginary part of the spectral transfer function is
calculated:
(
)
(
)
( )
( )
( )
( )
v vv
G j P jQ
ωω ω
= +
(8)
where:
(
)
( )
( )
( )
( ) ( ) ( )
22 4 3 3
1 2 11 2 2
4 3 3 22
11 2 2 1
22 2 22 2
1 1212 2 112
33
cos 2 2 cos 2 cos
22
33
cos 2 2 cos 2 cos cos
22
cos 4 cos cos 2 cos
2
vv v
v
v v vv
v vv v
vv
v
Pj
vv
vv
v
v vv
ππ
ω ω ω π ξωω ξωω
ω
ππ
ω π ξωω ξωω ω ω π
π
ωω π ξξωωω π ωω π ξωωω
++ +
=
+ + ++
+ ++
( ) ( )
( )
( )
2 22 2 2
1212 2 112 21 2
2 22
212 12
2
22 4 3 3 22
1 2 11 2 2 1
2
12 1 2
4 cos cos 2 cos 2 cos
22
2 cos
2
33
sin 2 2 sin 2 sin sin
22
4
vv v v
v
v v vv
v
vv
vv
v
vv
vv
ππ
ξξωωω π ωω π ξωωω ξωωω
π
ξωωω ωω
ππ
ωω ω π ξωω ξωω ωω π
ξξ ωωω
+
++ + +
+
++ + ++
( ) ( )
2
22 2 2
2 112 21 2
sin sin 2 sin 2 sin
22
vv v
vv
vv
ππ
π ωω π ξωωω ξωωω
++ +
(9)
645
( )
( )
( ) ( )
( ) ( ) ( )
22 4 3 3
1 2 11 2 2
4 3 3 22
11 2 2 1
22 2 22 2
1 1212 2 112
33
sin 2 2 sin 2 sin
22
33
cos 2 2 cos 2 cos cos
22
sin 4 sin sin 2 sin
2
vv v
v
v v vv
v vv v
vv
v
Qj
vv
vv
v
v vv
ππ
ω ω ω π ξωω ξ ωω
ω
ππ
ω π ξωω ξ ωω ω ω π
π
ωω π ξξωωω π ωω π ξωωω
++ +
=
+ + ++
+ ++
( ) ( )
( ) ( )
( )
2 22 2 2
1212 2 112 21 2
2
21 2
2
22 4 3 3 22
1 2 11 2 2 1
2
12 1 2
4 cos cos 2 cos 2 cos
22
2 sin
2
33
sin 2 2 sin 2 sin sin
22
4 sin
vv v v
v
v v vv
v
vv
vv
v
vv
vv
v
ππ
ξξωωω π ωω π ξωωω ξωωω
π
ξωωω
ππ
ωω ω π ξωω ξωω ωω π
ξξ ωωω π
+
++ + +
++ + ++
+
( )
2
22 2 2
2 112 21 2
sin 2 sin 2 sin
22
vv v
vv
v
ππ
ωω π ξωωω ξωωω
++
(10)
Knowing the real and imaginary part of the
spectral transfer function of the transducer, one can
determine the equation describing the logarithmic
amplitude step:
( )
( )
( )
( )
( )
( )
22
20
v vv
L log P Q
ω ωω
= +
(11)
and the equation describing the logarithmic phase
step:
( )
( )
( )
( )
( )
( )
( ) ( ) ( )
( )
43 3
11 2 2
43 3
11 2 2
22 2 22
1 12 1 2 2
22
1
33
sin 2 2 sin 2 sin
22
33
cos 2 2 cos 2 cos
22
sin 4 sin sin
cos 4
v
v
v
vv v
vv v
v vv
v
Q
arctg
P
vv
v
arctg
vv
v
v vv
v
ω
ϕω
ω
ππ
ω π ξωω ξωω
ππ
ω π ξωω ξ ωω
ωωπξξωωωπωωπ
ωω π ξ
= =
++ +
−
++ +
+ ++
+
( ) ( )
2 22
12 1 2 2
22
112 21 2
2 2 22
112 212 12
cos cos
2 sin 2 sin
22
2 cos 2 cos
22
vv
vv
vv
vv
vv
vv
ξωωωπωωπ
ππ
ξωωω ξωωω
ππ
ξωωω ξωωω ωω
++
+
++
(12)
Using the program written in the MATLAB
environment, which was used for conducting the
simulations of the equations describing the Bode plot
(11) and (12) of the membrane pressure transducer, a
response in the form of plots of logarithmic
magnitude and phase of the analyzed pneumatic
cascade was obtained. The plots are presented in Fig.
3 and Fig. 4.
The determined frequency response (Fig. 3 and
Fig. 4) correctly reflect the dynamics of the model. For
the parameter
1v =
, the logarithmic amplitude
response (Fig. 3) and phase response (Fig. 4) coincide
with the known responses of the 4th order oscillation
units. From the amplitude response (Fig. 3), one can
read the gain decrease, which is
80 /dB dek−
, and
from the phase response (Fig. 4), the phase shift
2
ϕπ
= −
for the parameter
1v =
, as it is in the classic
oscillation section of the 4th order.
The analysis of frequency responses (Fig. 3 and
Fig. 4) shows that the resonant pulsation depends on
the parameter
v
, and hence on the order of the
differential, in the differential equation describing the
studied system. By reducing the order, the resonant
pulsation increases. Hence, the smaller the phase shift
of the system is, the smaller the order of the
differential.
Figure 3. Logarithmic amplitude response of a pneumatic
cascade described by means of differential equation with
fractional derivatives of non-integer orders for the
parameter v in the range (0.8-1.2) [authors’ own elaboration]
Figure 4. Logarithmic phase response of a pneumatic
cascade described by means of differential equation with
fractional derivatives of non-integer orders for the
parameter v in the range (0.8-1.2) [authors’ own elaboration]
3 MATHEMATICAL MODEL OF A PNEUMATIC
MEMBRANE ACTUATOR DESCRIBED WITH A
FRACTIONAL CALCULUS
Figure 5. Pneumatic membrane actuator [11]
(Membrana - Membrane, Talerz oporowy - Retaining plate,
Sprężyna - Spring)
Based on Newton's laws, we can write:
( ) ( )
( ) ( )
2
0
2
dyt dyt
m R Cy t Ap t
dt
dt
+ +=
(13)
646
where:
[ ]
0
p Pa
- inlet pressure (forcing);
2
Am
- active surface of the membrane;
[ ]
m kg
- mass of the mobile system;
1
C Nm
−
- spring stiffness coefficient;
1
R Nsm
−
- viscous friction coefficient;
[
]
ym
- displacement of the actuator stem.
If the damping of the system
1
ξ
<
, then
2R mC<
.
Generalizing the equation (13), we get:
( ) ( )
( )
( )
2
00 0
RL v RL v
tt
m D yt R Dyt Cyt Ap t+ +=
(14)
where:
0v >
.
Applying the Laplace transform, assuming zero
initial conditions, for a fractional derivative of non-
integer order defined by Riemann - Liouville, we
obtain:
(
)
(
)
2
0
vv
RC
mY s s s Ap s
mm
+ +=
(15)
Thus the transfer function of non-integer order of
the membrane pneumatic actuator is obtained:
( )
( )
( )
( )
2
0
v
vv
A
Ys
m
Gs
RC
Ps
ss
mm
= =
++
(16)
In order to conduct a simulation, the data of a
pneumatic membrane actuator of the GMP brake T16
were used:
A - effective membrane surface: for r = 50mm
(diameter = 100mm) A = 0.00785 m2;
m - mass of the mobile system (membrane and
plunger); m = 0.12 + 0.2 = 0.32 kg;
C - spring stiffness coefficient; C = 1000N / m
R - viscous friction coefficient (resistance to
movement of moving parts) - R = 0,5Ns / m.
By substituting the above data for the equation
(16), the transfer function of the analyzed pneumatic
actuator was obtained:
( )
( )
( )
( )
2
0
0,025
1,56 3125
v
vv
Ys
Gs
Ps
ss
= =
++
(17)
The quantifier of the non-integer transfer function
has two complex roots. Therefore, you get:
( )
( )
( )
( )
( )
( )
0
1
0
0,025 3125 , 1,56; ,2
1
!
k
v
k
k
v
k
gt
t v v vk
k
ht
ε
∞
=
−
= − ++
∑
(18)
For the equation (8), you get the equations
describing the impulse and step response of the
incomplete order of the tested pneumatic actuator in
the form of:
( )
( )
( )
( )
( )
( )
( ) ( )
( )
( ) ( )
( )
2 11
,2
21
0
,2
1, 56
1
0,025 3125
!
1, 56
vk k
v
k
v
v v vk
k
v
vk k
v
k
v v vk
tE t
gt
k
ht
tE t
+−
∞
+
+
=
+
−
−
=
−
∑
(19)
whereas the function
( )
( )
,
k
Ez
αβ
is a special type of
Mittag-Leffler function.
By conducting a simulation of the pneumatic
transducer, a pulse and unit jump signal was applied
thus receiving the impulse and step response shown
in Figure 6 and Figure 7.
Figure 6. Impulse response of a pneumatic actuator
described with integer and non-integer order:
0,5
F
- for
0,5v
=
,
0,7
F
- for
0,7v =
,
0,9
F
- for
0,9v =
,
1,0
F
- for
1v =
,
2
C
- classic model (integer order) [authors' own
elaboration]
Figure 7. Step response of a pneumatic actuator described
with integer and non-integer order:
0,5
F
- for
0,5v =
,
0,7
F
- for
0,7v =
,
0,9
F
- for
0,9v =
,
1,0
F
- for
1v =
,
2
C
- classic model (integer order) [authors' own elaboration]
Fig. 6 and Fig. 7 show the impulse and step
response described with the formula (19) for selected
values
in the interval [0,1]. The impulse and step
responses of the tested pneumatic actuator described
by means of integer and non-integer differential
equations for the parameter
in the given scale
overlap. This indicates the correctness of the analyzed
actuator model. Fig. 6 and Fig. 7 show that for
increasing values of the orders of the analyzed
647
pneumatic actuator, impulse and step responses
behave like the second order oscillatory element. For
small orders that converge to 1, the answers behave
like the inertial element of the first order.
For the dependence (16) the spectral transfer
function of the actuator is obtained:
( )
( )
( )
( )
( )
( )
( ) ( )
2
2
cos sin cos sin
22
v
vv
v
vv
A
m
Gj
RC
jj
mm
A
m
Gj
RC
v jv v jv
mm
ω
ωω
ω
ππ
ω π πω
=
++
=
++ + +
(20)
By making elementary transformations, the real
and imaginary part of the spectral transfer function is
calculated, where:
(21)
Knowing the real and imaginary part of the
spectral transfer function of the transducer, we can
determine the equation describing the logarithmic
amplitude step:
( )
(
)
(
)
( )
( )
( )
22
20log
v vv
L PQ
ω ωω
= +
(22)
and the equation describing the logarithmic phase
step:
(
)
( )
( )
( )
( )
( )
( )
( )
2
2
sin sin
2
cos cos
2
vv
v
v
v
vv
Rv
v
Q
m
arctg arctg
R vC
P
v
mm
π
ω πω
ω
ϕω
π
ω
ω πω
+
= = −
++
(23)
Using the written program in the MATLAB
environment, simulations of equations (22) and (23)
were performed to obtain the frequency logarithmic
amplitude and phase responses of the actuator, which
are shown in Figure 8 and Figure 9.
Figure 8. Logarithmic amplitude response of a membrane
pneumatic actuator described by means of a differential
equation of non-integer orders for different sizes of the
parameter v (Siłownik - Pneumatic cylinder)
[authors' own elaboration]
Figure 9. Logarithmic phase response of a pneumatic
membrane actuator described by means of a differential
equation of non-integer orders for different sizes of the
parameter v (Siłownik - Pneumatic cylinder)
[authors' own elaboration]
The logarithmic frequency responses (Fig. 8 and
Fig. 9) show that for the parameter
1v =
, above the
resonant pulsation, the slope of the amplitude
response is
20 /dB dek−
, as it is in the second order
oscillating element. By reducing the order, the gain
decreases and the system behaves like the inertial
element of the first order.
The course of the logarithmic phase response
(Figure 9) confirms this trend. For the parameter
1v =
, the phase response overlaps the logarithmic
phase response of the classical second order
oscillating element (for the pulsation greater than the
resonant one, phase shift reaches the value
ϕπ
= −
).
By reducing the value of the order of the differential,
the system becomes the inertial element of the first
order, because for the pulsation greater than the
resonant one, phase shift decreases. For the
parameter
0,5v
=
, the phase shift will reach the
value
2
π
ϕ
= −
, as it is in the case of the inertial
element of the first order.
648
4 SUMMARY
The obtained responses, which arose from the
simulation of the dependencies resulting from the
solution of differential equations of integer orders,
overlap the responses of non-integer orders obtained
from the solution of differential equations of non-
integer order for the parameter
1v =
. This is
confirmed by the fact that the classical differential
calculus is a special case of the differential calculus of
any arbitrary order, and thus it proves that
mathematical models have been properly developed.
The use of the description of dynamic properties of
pneumatic systems based on the fractional calculus
will allow the authors of the article to analyze the
properties of a wide class of pneumatic systems of
arbitrary orders in the future.
BIBLIOGRAPHY
[1] Busłowicz M., Wybrane zagadnienia z zakresu liniowych
ciągłych układów niecałkowitego rzędu, Pomiary
Automatyka Robotyka nr 2/2010. (in Polish)
[2] Busłowicz M., Nartowicz T., Projektowanie regulatora
ułamkowego rzędu dla określonej klasy obiektów
z opóźnieniem, Pomiary Automatyka Robotyka, nr 2, s.
398-405, 2009. (in Polish)
[3] Kaczorek T., Wybrane zagadnienia teorii układów
niecałkowitego rzędu, Oficyna Wydawnicza Politechniki
Białostockiej, stron 271, ISSN 0867-096X, Białystok 2009.
(in Polish)
[4] Chwaleba A., Luft M., Właściwości i projektowanie
wybranych przetworników mechano–elektrycznych,
Zakład Poligraficzny Politechniki Radomskiej, Wyd. II
popr. i uzup., ISBN 83-88001-00-0, Radom 1998. (in
Polish)
[5] Luft M., Nowocień A., Cioć R., Pietruszczak D.,
Charakterystyki częstotliwościowe modelu
przetwornika ciśnienia opisanego równaniem
różniczkowym niecałkowitego rzędu, Logistyka nr
3/2015, ISSN 1231-5478, Poznań 2015. (in Polish)
[6] Luft M., Nowocień A., Pietruszczak D., Analiza
właściwości dynamicznych wybranych układów
pneumatycznych za pomocą rachunku różniczkowego
niecałkowitych rzędów. Część 1. Badania symulacyjne,
AUTOBUSY - Technika, Eksploatacja, Systemy
Transportowe; Eksploatacja i testy; ISSN 1509-5878, e-
ISSN 2450-7725, str. 1050-1055, Instytut Naukowo-
Wydawniczy "SPATIUM", AUTOBUSY 12(2017), Radom
2017. (in Polish)
[7] Luft M., Nowocień A., Pietruszczak D., Właściwości
dynamiczne wybranych podstawowych członów
automatyki niecałkowitych rzędów, AUTOBUSY -
Technika, Eksploatacja, Systemy Transportowe;
Eksploatacja i testy; ISSN 1509-5878, e-ISSN 2450-7725,
str. 1056-1060, Instytut Naukowo-Wydawniczy
"SPATIUM", AUTOBUSY 12(2018), Radom 2018. (in
Polish)
[8] Luft M., Nowocień A., Pietruszczak D., Analiza
właściwości dynamicznych wybranych układów
pneumatycznych za pomocą rachunku różniczkowego
niecałkowitych rzędów. Część 2. Badania laboratoryjne,
AUTOBUSY - Technika, Eksploatacja, Systemy
Transportowe; Eksploatacja i testy; ISSN 1509-5878, e-
ISSN 2450-7725, str. 1056-1060, Instytut Naukowo-
Wydawniczy "SPATIUM", AUTOBUSY 12(2017), Radom
2017. (in Polish)
[9] Luft M., Pietruszczak D., Nowocień A., Frequency
response of the pressure transducer model described by
the fractional order differential equation, TTS 12 (2016),
ISSN 1232-3829, Radom 2016.
[10] Luft M., Szychta E., Nowocień A., Pietruszczak D.,
Zastosowanie rachunku różniczkowo – całkowego
niecałkowitych rzędów w matematycznym
modelowaniu przetwornika ciśnienia, Autobusy
nr 6/2016, ISSN 1509-5878, Instytut Naukowo-
Wydawniczy SPATIUM, Radom 2016 (in Polish)
[11] Nowocień A., Analiza właściwości dynamicznych
układów pneumatycznych za pomocą rachunku
różniczkowego niecałkowitych rzędów, Rozprawa
doktorska, Biblioteka Główna Uniwersytetu
Technologiczno- Humanistycznego im. Kazimierza
Pułaskiego w Radomiu, Radom 2017, (Promotor: Prof.
dr hab. Inż. Mirosław Luft; Promotor pomocniczy: dr
inż. Daniel Pietruszczak (in Polish)
[12] Nowocień A., Luft M., Pietruszczak D., Zastosowanie
rachunku różniczkowo – całkowego niecałkowitych
rzędów w nauce i technice. Logistyka nr 3/2014. (in
Polish)
[13] Ostalczyk P., Zarys rachunku różniczkowo-całkowego
ułamkowych rzędów. Teoria i zastosowanie w
automatyce, Wydawnictwo Politechniki Łódzkiej, stron
430, ISBN 978-83-7283-245-0, Łódź 2008. (in Polish)
[14] Pietruszczak D., Analiza właściwości układów
pomiarowych wielkości dynamicznych z
wykorzystaniem rachunku różniczkowo – całkowego
ułamkowych rzędów, Rozprawa doktorska,
Biblioteka Główna Uniwersytetu Technologiczno-
Humanistycznego im. Kazimierza Pułaskiego w
Radomiu, Radom 2013. (in Polish)
[15] Podlubny I., Fractional Differential Equations. An
Introduction to Fractional Derivatives, Fractional
Differential Equations, Some Methods of Their Solution
and Some of Their Applications, Academic Press,
368 pages, ISBN 0125588402, San Diego-Boston-
New York-London-Tokyo-Toronto 1999.
[16] Mościński J., Ogonowski Z. (red.), Advanced Control
with MATLAB and SIMULINK, Pearson Higher
Education, 272 pages, ISBN 013309667X, 1995.
[17] Rudra P., (tłumacz: Korbecki M.), Matalb 7 dla
naukowców i inżynierów, Wydawnictwo Naukowe
PWN, stron 280, ISBN 9788301160579, Warszawa 2010.
[18] http://www.mathworks.com – strona internetowa
producenta programu MATLAB
