MODELOWANIE LOTU SAMOLOTU TRANSPORTOWEGO Z UWZGLĘDNIENIEM ZMIENNYCH OBCIĄśEŃ ATMOSFERYCZNYCH

MODELOWANIE INśYNIERSKIE ISSN 1896-771X 36, s. 193-200, Gliwice 2008

MODELOWANIE LOTU SAMOLOTU TRANSPORTOWEGO Z UWZGLĘDNIENIEM ZMIENNYCH OBCIĄśEŃ

ATMOSFERYCZNYCH

G

RZEGORZ

K

OWALECZKO

, J

AROSŁAW

K

RZONKALLA

, M

IROSŁAW

N

OWAKOWSKI

Instytut Techniczny Wojsk Lotniczych e-mail:g.kowaleczko@chello.pl

Streszczenie. W pracy przedstawiono sposób modelowania lotu samolotu transportowego w burzliwej atmosferze. Problem rozwiązano, wykorzystując układ równań róŜniczkowych opisujących ruch przestrzenny samolotu.

Turbulencje atmosferyczne określono z uwzględnieniem ich stochastycznego charakteru. Siły i momenty aerodynamiczne generowane przez skrzydło wyznaczono poprzez numeryczne całkowanie odpowiednich wyraŜeń. Określono reakcję samolotu na zaburzenia pochodzące od turbulencji.

1. WSTĘP

W typowych analizach z zakresu dynamiki lotu statków powietrznych zazwyczaj przyjmuje się, Ŝe lot odbywa się w spokojnej, nieruchomej atmosferze - wpływ ruchu mas powietrza jest pomijany [1-4]. W rzeczywistości porywy wiatru mogą w zasadniczy sposób wpływać na realizację zadań lotnych bezpośrednio wpływając równieŜ na bezpieczeństwo załogi.

W pracach dotyczących oddziaływania wiatru na statki powietrzne często stosuje się modelowanie pola wiatru w postaci zaleŜności funkcyjnych przybliŜających rzeczywiste zjawiska atmosferyczne jedynie z ograniczoną dokładnością [5] - wiadomym jest, Ŝe mają one charakter stochastyczny, co w tym przypadku jest pomijane. Uwzględnienie takiego charakteru wiatru wymaga odpowiedniej modyfikacji opisu ruchu samolotu. Dotyczy to w szczególności procedur wyznaczania sił i momentów aerodynamicznych. JeŜeli rozmiary samolotu są niewielkie w porównaniu z wymiarem charakterystycznym oddziałujących na niego podmuchów, to moŜna załoŜyć, Ŝe wiatr zmienia kąt natarcia i kąt ślizgu wszystkich powierzchni nośnych samolotu o tę samą wartość. Zmienione kąty moŜna wyznaczyć, znając wektor prędkości samolotu względem powietrza. Prędkość ta jest równa róŜnicy prędkości samolotu względem Ziemi i prędkości wiatru względem Ziemi.

Omówione powyŜej przybliŜenie jest słuszne dla samolotów małych, w szczególności dla bezpilotowych statków powietrznych [6]. JeŜeli wymiary samolotu są duŜe, to turbulencje mogą być róŜne w rejonie róŜnych części samolotu. W związku z tym dokonano modyfikacji sposobu opisu dynamiki ruchu samolotu. Polega ona na uwzględnieniu zmienności podmuchu wiatru wzdłuŜ rozpiętości skrzydła, co powoduje powstanie róŜnych kątów natarcia w kolejnych przekrojach skrzydła. Zastosowanie takiego modelowania oznacza konieczność

wyznaczania obciąŜeń aerodynamicznych poprzez całkowanie odpowiednich wyraŜeń wzdłuŜ skrzydła w kolejnych chwilach czasu.

Analizując wpływ burzliwości atmosfery na lot samolotu, kluczowym zagadnieniem, które naleŜy rozwiązać, jest odtworzenie stochastycznej struktury pola wiatru. W prowadzonych badaniach wykorzystano model Shinozukiego [7], w którym moc widmową podmuchu określono stosując spektrum Drydena.

Zaproponowana metoda pozwala na wyznaczenie przebiegów parametrów lotu w funkcji czasu. Daje ona teŜ moŜliwość oceny zmienności obciąŜeń aerodynamicznych wzdłuŜ skrzydła w wybranych fazach lotu.

2. MODEL RUCHU SAMOLOTU

ZałoŜono, Ŝe samolot jest bryłą sztywną o sześciu stopniach swobody, której masa i momenty bezwładności są niezmienne. Ma on płaszczyznę symetrii geometrycznej, masowej i aerodynamicznej.

2.1. Układy współrzędnych i ich transformacja

Dla tak zdefiniowanego modelu określone zostały układy współrzędnych, które wykorzystano do opisu ruchu samolotu. Są to: 1. układ związany z Ziemią o początku w środku masy samolotu Oxgygzg, 2. układ związany z samolotem Oxyz o początku w środku masy samolotu, 3. układ związany z przepływem Oxayaza o początku w środku masy samolotu. Układy te powiązane są ze sobą poprzez kąty przejścia pozwalające, poprzez kolejne obroty, na pokrycie się jednego układy z drugim. Pokazano je na rys.1 i 2.

xg

yg

zg

x

y

z

Rzut osi Ox na płaszcznę Ox ygg

x’

y’

z”

xa

ya

za=z’

x

y=y’

z

Rzut osi Ox na płaszcznę Ox ygg

x’

Rys.1 Układy współrzędnych Oxyz i Ox_gy_gz_g oraz kąty przejścia pomiędzy nimi

Rys.2. Układy współrzędnych Oxyz i Ox_ay_az_a oraz kąty przejścia pomiędzy nimi

Macierze transformacji są odpowiednio równe:

Macierz transformacji L_{s /g} z układu O_gx_gy_gz_g do układu Oxyz:

















Φ Θ Φ

Ψ

− Φ Θ Ψ Φ

Ψ + Φ Θ Ψ

Φ Θ Φ

Ψ + Φ Θ Ψ Φ

Ψ

− Φ Θ Ψ

Θ

− Θ

Ψ Θ

Ψ

=

cos cos sin

cos cos

sin sin sin

sin cos

sin cos

sin cos cos

cos sin

sin sin cos

sin sin

sin cos

sin cos

sin cos

cos

/ g

Ls (1)

Macierz transformacji L_s/a z układu Ox_ay_az_a do układu Oxyz :

















−

−

−

=

α β

α β

α

β β

α β

α β

α

cos sin

sin cos

sin

0 cos

sin

sin sin

cos cos

cos

/ a

Ls (2)

Przeliczenia z układu Oxyz do układu O_gx_gy_gz_g oraz do układu Ox_ay_az_a realizowane są poprzez macierze:

1 / /

= −_{s g} s

g L

L , L_a/s =L⁻_s¹_/a (3)

2.2. Równania ruchu

Równania ruchu postępowego i równania ruchu obrotowego mają następującą postać wektorową:

V =F dt

md , _gir

dt

dK M M

+

= (4) gdzie: m – masa samolotu; V - prędkość środka masy samolotu w inercjalnym układzie współrzędnych; F - wypadkowy wektor sił działających na samolot; K - wektor momentu pędu (kręt) samolotu; M - wypadkowy moment sił działających na samolot; Mgir - moment giroskopowy.

Przyjmując, Ŝe w układzie Oxyz poszczególne wektory mają następujące składowe:

V=[u,v,w]^T , F=

[

^F_x^,F_y^,F_z

]

^T, K=

[

^K_x^,K_y^,K_z

]

^T,

[

^L,M,N

]

^T

=

M , Mgir =

[

Lgir,Mgir,Ngir

]

^T, Ω=[p,q,r]^T

(5) równania (5) i (6) moŜna zapisać w postaci:

z y x

F qu pv w m

F pw ru v m

F rv qw u m

=

− +

=

− +

=

− +

) (

) (

) (

&

&

&

( ) ( )

( ) ^{( )}

( ) (

_{x y}

)

_gir

xz z

gir x

z xz

y

gir z

y xz

x

N N pq I I qr p I r I

M M rp I I p r I q I

L L qr I I pq r I p I

+

=

−

−

−

−

+

=

−

−

−

−

+

=

−

− +

−

&

&

&

&

&

2 2

(6)

Uzupełnia się je związkami kinematycznymi dotyczącymi kątów oraz współrzędnych środka masy samolotu:

Φ& = p+

(

rcosΦ+qsinΦ

)

tgΘ

Θ& =qcosΦ−rsinΦ

(

^{Φ+ Φ}

)

= Θ

Ψ cos sin

cos

1 r q

&

















=

















=

















w v u L w v u

z y x

s g g g g

g g g

/

&

&

&

(7)

W efekcie otrzymuje się zamknięty układ równań róŜniczkowych zwyczajnych, który pozwala obliczyć przebiegi parametrów lotu X=

[

u,v,w,p,q,r,Φ,Θ,Ψ,xg,yg,zg

]

^T.

2.3. Siły działające na samolot

Siła F równa jest sumie sił: - aerodynamicznej R, - cięŜaru samolotu Q, - ciągu zespołu napędowego P.

P Q R

F= + + (8) Siła aerodynamiczna R jest równa sumie sił powstających na elementach konstrukcyjnych samolotu:

V H sk

k R R R

R

R = + + + (9) Poszczególne indeksy oznaczają odpowiednio: k – kadłub, sk – skrzydło, H – usterzenie poziome, V – usterzenie pionowe.

Składowe tej siły określane są zwykle w układzie związanym z przepływem Ox_ay_az_a: V S

C P

R_{x xa xa}

a 2

2

ρ *

−

=

−

= , V S

C P

R_{y ya ya}

a 2

2

ρ *

=

= , V S

C P

R_{z za za}

a 2

2

ρ *

−

=

−

= (10)

Pxa jest siłą oporu, Pza – siłą nośną, Pya – siłą boczną. V_* – prędkość samolotu względem powietrza. CięŜar samolotu Q ma w układzie Oxgygzg jedną składową wzdłuŜ osi Ozg:

[

0,0,mg

]

^T

=

Q . Ciąg zespołu napędowego P ma kierunek i zwrot zgodny z osią Ox samolotu.

Składowe siły F w układzie związanym z samolotem są odpowiednio równe:















 +















 +

















=

















0 0 0

0

/ /

P

mg R

R R

F F F

g s z

y x a s z y x

a a a

L

L (11)

Podobnie jak siła aerodynamiczna, moment aerodynamiczny M powstaje w wyniku opływu przez powietrze kadłuba, skrzydła oraz usterzenia poziomego i pionowego. Przyjęto, Ŝe jest ona równy sumie momentów od sił powstających na tych elementach konstrukcyjnych:

V H sk

k M M M

M

M= + + + (12) Składowe tego momentu w układzie Oxyz związanym z samolotem są następujące:

V Sl C

L _l

2

2

ρ *

= , _m V Sb_a

C

M 2

2

ρ *

= , V Sl

C

N _n

2

2

ρ *

= (13) L jest momentem przechylającym, M – momentem pochylającym, N – momentem odchylającym, l – rozpiętość skrzydła, ba – średnia cięciwa aerodynamiczna skrzydła.

Moment giroskopowy M_gir powstaje, jeŜeli na samolocie wykonującym ruchy obrotowe zabudowane są wirujące masy. Jest on równy:

∑

=

×

= ^m

k

k k gir

1

Ω ω J

M (14) gdzie: J – moment bezwładności k-tej wirującej masy;_k ω – wektor prędkości obrotowej k-_k tej wirującej masy, m – ilość wirujących mas, które uwzględniono.

2.4. Siły i momenty aerodynamiczne skrzydła

Jak napisano powyŜej, siły i momenty aerodynamiczne skrzydła oblicza się poprzez obliczenie odpowiednich całek wzdłuŜ skrzydła samolotu.

∫

= _x

sk

x dR

R , ^{Rysk =}

∫

^{dRy, Rzsk =}

∫

^dRz (15)

∫

= _sk

sk L

L , ^{Msk =}

∫

^{dMsk , Nsk =}

∫

^dNsk (16)

WyraŜenia podcałkowe oblicza się na podstawie znajomości elementarnych sił powstających w kolejnych przekrojach skrzydła:

ρ κ ρ α

ρ α V b y d

C V dS

C dP

dR _{xa xapr P} ^P _{xapr P} ^P _P

a ( )

) 2 2 (

) (

2

* 2

* =−

−

=

−

= (17)

ρ κ ρ α

τ α V b y d

C V dS

C dP

dR _{za zapr P} ^P _{zapr P} ^P _P

a ( )

) 2 2 (

) (

2

* 2

* =−

−

=

−

= (18)

P z

sk dR y

dL =− ⋅ , dM_sk =dR_z⋅x_P −dR_x⋅z_P, dN_sk =dR_x⋅y_P (19) xP oraz y_P są współrzędnymi przekroju skrzydła w układzie związanym z samolotem.

W obliczeniach stosuje się transformację:

















=

















=

















a a

dR dR dR

dR dR dR

dR dR

pa p p s p

s z y x

τ ρ

τ κ ρ

/ 0

/

/ L L

L (20)

linia 1/4 cię

ciw rzut lin

ii 1/4 cięciw na płaszczyznę Oxyz

y_P x

P

zP

x ψ y χ

z P

O’

κ

τ ρ

z y

x O

,

,

,

P*

V+VΩ

Vw

V

w

κ κ v

τ

α ρ

u β

P

P*

P* P*

P

P

ρa

τa a

αP

Rys.3. PołoŜenie układu związanego z profilem skrzydła względem układu Oxyz

Rys.4. Określenie kąta natarcia profilu αP i kąta „ślizgu” βP

Macierze transformacji są równe:

















−

=

ψ ψ

ψ ψ

cos sin

0

sin cos

0

0 0

1

/ s

Lp ,















 −

=

P P

P P

pa p

α α

α α

cos 0 sin

0 1 0

sin 0

cos

L / (21)

Występujące tu kąty pokazano na rys. 3 i 4. Kąt natarcia przekroju αP zaleŜy od rozkładu prędkości w tym przekroju:

*

arctan * P

P

P u

= w

α

(22)

*

wP i w_P* są składowymi prędkości przekroju względem powietrza V_P*. Prędkość ta jest równa róŜnicy pomiędzy prędkością bezwzględną punktu P V_P i prędkością wiatru Vw:

w w

P

P V V V V V

V _* = − =( + _Ω)− (23) Prędkość V_P jest sumą prędkości środka masy samolotu V i prędkości wynikającej z ruchu obrotowego V_Ω.

3. MODEL POLA WIATRU

W celu modelowania losowych przebiegów czasowych turbulencji wiatru zastosowano metodykę opisaną szerzej w pracach [60, 61]. ZałoŜono tam, Ŝe proces losowy moŜe być opisany jako superpozycja funkcji okresowych, których częstość oscylacji zmienia się w sposób losowy, zaś amplituda zaleŜy od gęstości widmowej mocy.

Podstawowe wyraŜenie pozwalające obliczyć pole wiatru jako proces stochastyczny ma następującą postać:

( ) (

l jl

)

i

j L

l

l ij

i H φ

υ =

∑∑

∆ +

= =

r Ω Ω Ω

r ^'

1 1

cos 2

)

( (24)

Ω - wektor częstości „przestrzennej”; r – wektor określający połoŜenie rozpatrywanego punktu, φ_jl - wzajemnie niezaleŜne, losowo zmienne przesunięcia fazowe o wartościach z przedziału 0÷2π. H jest amplitudą związaną z macierzą gęstością widmową mocy ΦΦΦ Φ zaleŜnością:

) ( ) ( )

(Ω H Ω H Ω

Φ = ⋅ ^T (25) Macierz gęstości widmowej moŜna określić, wykorzystując jeden z dostępnych w literaturze modeli. Tu zastosowano model Drydena:

( )

[ ]

( )

( )

( )

^















Ω + Ω Ω

+ Ω + Ω

Ω

−

Ω Ω

− Ω

+ Ω +

Ω + Ω +

= Ω Ω

2 2 2 2 2 2 2

2 2

2 2

2 5 2 2 2 2 2

3 0

0

0 4

1 3

0 3

4 1

4 1 ) , (

y x w y x w w

y x

w y x y

x w

y x w w y x

L L

L

L L

L

L σ

Φ π (26)

W przeprowadzonych analizach uwzględniono hipotezę Taylora, zgodnie z która pole wiatru, w którym porusza się samolot moŜe być traktowane jako niezmienne w czasie. Zakładając dodatkowo, Ŝe nie zaleŜy ono od wysokości lotu, a jedynie od współrzędnych x_g i y_g, składowe wiatru są równe:

( ) [ ]

∑∑

= =

+ Ω + Ω

∆Ω

⋅

∆Ω

⋅ Ω Ω

= ^x

x y

y

y x y

x y

x

L

l L

l

l l g yl g xl y

x yl

xl g

g

wg x y H x y

u

1 1

1 ,

, 11

t( , ) , 2 cos φ (27)

( ) [ ]

( ) [ ]

∑∑

∑∑

= =

= =

+ Ω + Ω

∆Ω

⋅

∆Ω

⋅ Ω Ω +

+ + Ω + Ω

∆Ω

⋅

∆Ω

⋅ Ω Ω

=

x

x y

y

y x y

x y

x x

x y

y

y x y

x y

x

L

l L

l

l l g yl g xl y

x yl

xl L

l L

l

l l g yl g xl y

x yl

xl g

g wg

y x

H

y x

H y

x v

1 1

2 ,

, 22

1 1

1 ,

, 21

t

cos 2

,

cos 2

, )

, (

φ φ

(28)

( ) [ ]

( ) [ ]

( ) [ ]

∑∑

∑∑

∑∑

= =

= =

= =

+ Ω + Ω

∆Ω

⋅

∆Ω

⋅ Ω Ω +

+ + Ω + Ω

∆Ω

⋅

∆Ω

⋅ Ω Ω +

+ + Ω + Ω

∆Ω

⋅

∆Ω

⋅ Ω Ω

=

x

x y

y

y x y

x y

x x

x y

y

y x y

x y

x x

x y

y

y x y

x y

x

L

l L

l

l l g yl g xl y

x yl

xl L

l L

l

l l g yl g xl y

x yl

xl L

l L

l

l l g yl g xl y

x yl

xl g

g wg

y x

H

y x

H

y x

H y

x w

1 1

3 ,

, 33

1 1

2 ,

, 32

1 1

1 ,

, 31

t

cos 2

,

cos 2

,

cos 2

, )

, (

φ φ

φ

(29)

4. PRZYKŁADOWE WYNIKI MODELOWANIA

W wyniku zastosowania omówionego powyŜej sposobu modelowania burzliwości otrzymano szereg róŜnych pól pola wiatru, które spełniają wymóg niepowtarzalności i mają charakter stochastyczny. Przykładowe rozkłady dwóch składowych turbulentnych wiatru pokazano na rys. 5. Przeprowadzone obliczenia sprawdzające wykazały, Ŝe w kaŜdej kolejnej symulacji otrzymywano róŜne rozkłady.

Rys.5. Rozkład składowych turbulencji wiatru w płaszczyźnie poziomej

W czasie lotu samolotu wszystkie parametry lotu ulegały zaburzeniom. Zaobserwowano dwa typy przebiegów. Pierwszy dotyczył ruchów określanych w literaturze jako fugoidalne (długookresowe). Reprezentują je pokazane na rys.6 i 7 prędkość lotu i kąt pochylenia samolotu. Taka reakcja samolotu jest odmienna od badanych wcześniej zachowań małego samolotu bezpilotowego BSL, którego prędkość lotu zmieniała się współbieŜnie z turbulencją atmosfery. Natomiast na rys.8 pokazano przykładowe zmiany kąta natarcia samolotu, którego dotyczą ruchy krótkookresowe. Widać, Ŝe w tym przypadku przebieg ten jest zgodny ze zmianami prędkości pionowej określonymi w środku masy samolotu. Podobne wyniki otrzymano równieŜ dla BSL.

60 65 70 75 80 85

0 50 100 150 200 250 300 350 400 450 500

t [s]

V [m/s]

-15 -10 -5 0 5 10 15

0 50 100 150 200 250 300 350 400 450 500

t [s]

ΘΘΘΘ [deg]

Rys.6. Przebieg prędkości lotu samolotu Rys.7. Przebieg kąta pochylenia samolotu

-2 -1 0 1 2 3 4 5 6

0 50 100 150 200 250 300 350 400 450 500

t [s]

αααα [deg]

-5 -4 -3 -2 -1 0 1 2 3 4 5

0 50 100 150 200 250 300 350 400 450 500

t [s]

Vind_z [m/s]

Rys.8. Przebieg kąta natarcia w środku masy Rys.9. Przebieg składowej pionowej wiatru

2.030 2.035 2.040 2.045 2.050 2.055 2.060 2.065 2.070 2.075 2.080

-12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12

y [m]

α [deg]

3375 3380 3385 3390 3395 3400 3405 3410 3415 3420

-12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12

y [m]

Fττττ [N/m]

Rys.10. Rozkład kąta natarcia wzdłuŜ rozpiętości skrzydła

Rys.11. Rozkład obciąŜenia normalnego wzdłuŜ rozpiętości skrzydła

Zastosowany sposób modelowania pozwala wyznaczyć równieŜ rozkład kąta natarcia i obciąŜeń aerodynamicznych wzdłuŜ skrzydła. Pokazano to rys.9 i 10. Parametry te zmieniają się w kaŜdej chwili lotu. Analiza takich przebiegów umoŜliwia oszacowanie maksymalnych sił i momentów oraz widma obciąŜeń, które mogą oddziaływać na konstrukcję skrzydła. MoŜe to być waŜne dla oceny wytrzymałości zmęczeniowej skrzydła.

LITERATURA

1. Etkin B.: Dynamics of atmosphere flight. New York : Johson Willey & Sons Inc., 1972.

2. Kowaleczko G.: Zagadnienie odwrotne w dynamice lotu statków powietrznych.

Warszawa : Wyd. WAT, 2003.

3. Ostosławskij N. W.: Aerodynamika samoleta. Moskwa: Oborongizdat, 1957.

4. Stevens B., Lewis F.:Aircraft control and simulation. New York : Johson Willey & Sons Inc., 1992.

5. Holbit F.M.: Gust loads om aircraft: concepts and applications. “AIAA Education Series”

Washington D.C., 1988.

6. Czechowicz B., Hajduk J., Kowaleczko G., Nowakowski M.: Numeryczna analiza dynamicznych własności bezpilotowego statku powietrznego HOB-bit W; Materiały II Międzynarodowej Konferencji „Naukowe Aspekty Bezzałogowych Obiektów Latających”. Kielce 2006.

7. Shinozuka M.: Simulation of multivariate and multidimensional random processes.

“Journal of the Acoustical Society of America” 1971, Vol. 49.

MODELLING OF A FLIGHT OF A TRANSPORT AIRCRAFT INCLUDING CHANBEABLE ATMOSPHERIC LOADS

Summary. A method of modelling of a flight for a transport aircraft performing flight in a turbulent atmosphere is presented. This problem was resolved on the basis of the set of differential equations, which describes a spatial motion of the plane. Atmospheric turbulences were described taking into account their stochastic character. Aerodynamic forces and moments produced by a wing were determined by numerical integration of appropriare expressions. A reaction of flight to disturbances generated by turbulences was determined.