WPŁYW PRĘDKOŚCI LOTU ENTOMOPTERA NA GENEROWANE OBCIĄŻENIA
AERODYNAMICZNE
Paweł Czekałowski
1a, Krzysztof Sibilski
1b1Zakład Inżynierii Lotniczej, Politechnika Wrocławska
e-mail: apawel.czekalowski@pwr.wroc.pl, bkrzysztof.sibilski@pwr.wroc.pl
Streszczenie
Przedstawiony artykuł ściśle nawiązuje do prac nad stworzeniem obiektu latającego, który generuje siłę nośną tak jak latające owady, czyli entomoptera. Przedstawiona praca dotyczy badań doświadczalnych prowadzonych na me- chanizmie trzepoczącym pracującym w tunelu wodnym. Obiekt wyposażony jest w dwa skrzydła. Każde z nich może wykonywać dowolny ruch kulisty (obrót wokół punktu). Celem pracy było wyznaczenie maksymalnych pręd- kości lotu i uzyskiwanych wartości sił nośnych obiektu dla różnych sposobów ruchu skrzydeł. Eksperyment polegał na pomiarze sił hydrodynamicznych generowanych przez obiekt za pomocą wagi tensometrycznej. Ruch skrzydeł odbywał się w jednej płaszczyźnie. Parametrami doświadczenia były prędkość postępowa oraz kąt pochylenia płaszczyzny trzepotania. Uzyskane wyniki pozwalają na przewidzenie maksymalnej prędkości lotu rzeczywistego obiektu dla zadanego udźwigu, wskazują ewentualne sposoby poprawy osiągów. Dają również istotne informacje z punktu widzenia mechaniki lotu, a konkretnie stateczności podłużnej układu trzepoczącego. Artykuł oprócz ana- lizy wyników eksperymentu zawiera także opis metodologii pomiarów
INFLUENCE OF CRUISE FLIGHT SPEED OF ENTOMOPTER ON AERODYNAMICS LOADS
Summary
The article is close connected with building flying object, that fly like an insect (entomopter). Present work con- cerns on experimental investigations conducted on flapping mechanism in water tunnel. The object is dipterous.
Each wing can perform various spherical motions (wing is rotated around point). The aim of work is to find max- imal cruise flight speed and the way, how the lift force is changing due to this speed. The experiment consist of measurements of hydrodynamic forces generated by object. The motion of the wing in this case was two dimen- sional, it was rotated around two axis. Forward velocity and angle of inclination of flapping plane were parameters of the experiment. Test was performed for various pitch angle trajectories. In addition, pitching moment and lift force derivatives of the forward velocity were designated.
1. SPIS OZNACZEŃ
A,B - współczynnik szeregu Fouriera c -cięciwa
CL - współczynnik siły nośnej f - częstotliwość trzepotania
FIn - składowa siły bezwładności działająca w kierunku siły normalnej
FIt - składowa siły bezwładności działająca w kierunku siły stycznej
J - posuw względny
ms - masa skrzydła
N - składowa normalna siły hydrodynamicznej
V - siła prostopadła do wektora prędkości postępowej (nośna)
H - siła równoległa do wektora prędkości postępowej (napędowa)
T - składowa styczna siły hydrodynamicznej u - prędkość postępowa
PA - kąt pochylenia
Re - liczba Reynoldsa
R - odległość końcówki skrzydła od osi obrotu R0 - odległość nasady skrzydła od osi obrotu r - odległość sekcji skrzydła od osi obrotu
γ φ,
r
sc - odległość środka ciężkości skrzydła od odpo- wiednich osi obrotuγ φ,
I
- momenty bezwładności względem odpowiednich osiΦ - kątowy zakres ruchu względem średniej płaszczyzny γ -kąt nastawienia
φ
- położenie azymutalne skrzydłaλ - wydłużenie skrzydeł
S R2
= 4
λ (S-sumaryczna powierzchnia skrzydeł)
ρ -gęstość wody
υ - kinematyczny współczynnik lepkości ω - prędkość kątowa skrzydła
γ
ε
φ, - przyspieszenie kątowe skrzydła względem odpo- wiednich osi2. TŁO I GENEZA PROBLEMU
Przedstawiony eksperyment przeprowadzony został w ramach projektu rozwojowego: „System monitorowa- nia i detekcji aktywności istot żywych w pomieszcze- niach zakrytych i budynkach oraz ochrony i monitoro- wania obiektów z latającym nanorobotem w układzie entomoptera” finansowanego ze środków NCBiR. Pro- jekt ma na celu zbadanie możliwości budowy i wykorzy- stania obiektów klasy mikro w układzie entomoptera.
Eksperyment był doświadczalną próbą sprawdzenia w jaki sposób badany obiekt będzie zachowywał się podczas lotu postępowego. Doświadczenie przede wszystkim miało odpowiedzieć na pytania: jak szybko będzie w stanie polecieć oraz w jaki sposób należałoby nim sterować?
Modelowanym obiektem był prototypowy entomopter o konstrukcji przedstawionej na rys. 1. Podstawowe parametry charakteryzujące kinematykę ruchu zostały zestawione w tabeli 1. Obiekt jest dwuskrzydły. Do napędu służy miniaturowy silnik elektryczny, przekaza- nie mocy na skrzydła realizowane jest przez układ korbowy. Rozpiętość skrzydeł wynosi 200mm. Całkowi- ty zakres ruchu głównego to 140o. Częstotliwość trzepo- tania zawiera się w zakresie 10-20Hz, a średnia liczba
Reynoldsa charakteryzująca opływ skrzydeł to przedział 11000-23000.
Tabela 1. Parametry entomoptera
Rys.1. Projekt entomoptera
3. IDEA EKSPERYMENTU
Ideą eksperymentu było zmierzenie sił hydrodynamicz- nych generowanych przez skrzydła robota dla różnych prędkości postępowych i różnych kątów pochylenia obiektu (względem wektora prędkości) i określenie na tej podstawie maksymalnej prędkości lotu postępowego.
Eksperyment przeprowadzony został w tunelu wodnym.
Opis aparatury badawczej można znaleźć w [1] oraz [2]
Powodem, dla którego warto jest przeprowadzać tego typu testy w wodzie, jest możliwość radykalnego zmniejszenia częstotliwości trzepotania, a w efekcie zminimalizowanie udziału w wynikach sił i momentów bezwładności, które utrudniają (a mogą wręcz uniemoż- liwić) identyfikacje sił aerodynamicznych. Tego typu badania bardzo często prowadzone są właśnie w tym ośrodku, czego przykładem mogą być artykuły [3], [4], [5]. Jako model entomoptera do badań w tunelu wod- nym wykorzystany został specjalnie do tego zaprojekto- wany robot. Opis jego konstrukcji można znaleźć w pozycji [6].
Ze względów technicznych testy wykonywane były ze średnią liczbą Reynoldsa 7500 (trudności ze zwiększe- niem częstotliwości, ograniczenie rozpiętości skrzydeł ze względu na szerokość sekcji pomiarowej). Liczba ta stosowana jako kryterium podobieństwa dynamicznego opływu zdefiniowana została w następujący sposób:
masa M 10g
odległość końcówki skrzy-
dła od osi obrotu R 100mm
amplituda ruchu Φ 140o
częstotliwość trzepotania f 10-20Hz
liczba Reynoldsa Re 11000-23000
(1)
Mechanizm umożliwia obracanie każdym ze skrzydeł względem trzech osi, jednakże na potrzeby eksperymen- tu napęd jednej z osi został wyłączony. Zmiany kątowe- go położenia skrzydła podczas jednego cyklu przedsta- wione zostały na rys.6 i 7. Rys. 2. przedstawia definicję kątów opisujących położenie skrzydła. Przebieg zmian kąta azymutalnego odpowiada zaprojektowanemu układowi korbowemu, natomiast kąta nastawienia został przyjęty arbitralnie. Analitycznie ruch opisany jest za pomocą szeregu Fouriera:
[ ]
01
) 2
sin(
) 2
cos(
) ( ),
(t t A i f t B i f t A
N
i
i
i ⋅ ⋅ ⋅ + ⋅ ⋅ ⋅ +
=
∑
=
π π
γ φ
. (2)[7]
Współczynniki zostały zestawione w tabeli 2. Kształt obrysu skrzydła (rys. 3) został zaczerpnięty z przyrody (skrzydło muchy domowej). Przekrój stanowi płaska płytka zaostrzona na krawędzi natarcia i spływu.
Tabela 2. Współczynniki szeregu Fouriera
Rys. 3 Kształt skrzydła
Rys. 2. Definicja położenia skrzydła
Zakres sprawdzonych prędkości postępowych to u=0- 0,08604m/s. Prędkość postępową można wyrazić w sposób bezwymiarowy poprzez posuw względny, czyli stosunek prędkości postępowej do średniej prędkości końcówki skrzydła.
R f J u
⋅
⋅
= Φ
2
(3) [7]Wspomniana wielkość zapożyczona została z nazewnic- twa lotniczego, która służy do określania warunków pracy śmigła. W przypadku entomoptera ruchu skrzydła również jest okresowy, więc możliwy jest analogiczny sposób opisu warunków pracy. Przy posługiwaniu się tą bezwymiarową wielkością zakres wymienionych prędko- ści odpowiada wartościom J=0-0,54. Zakres prędkości ustalono na podstawie danych statystycznych (maksy- malny posuw względny osiągany przez latające owady [8]). Kąt pochylenia zdefiniowany jako kąt pomiędzy osią silników a kierunkiem prostopadłym do wektora prędkości postępowej, zmieniany był w zakresie 0-40o.
4. METODOLOGIA
Do pomiaru siły wykorzystana została pięcioskładniko- wa waga tensometryczna będąca na wyposażeniu tunelu.
Jedynie odczyty z sekcji odpowiadających za siłę nor- malną oraz moment pochylający były brane pod uwagę.
Pozostałe wskazują jedynie szum oraz siły i momenty będące efektem asymetrii ustawienia obiektu. Model wraz z wagą umieszczony został na suporcie tunelu wodnego, który umożliwił zmianę kąta pochylenia podczas cyklu badań. Lot z prędkością postępową imitowany był przez przepływ wody w tunelu. Każdy
φ
γi A B A B
0 0 90
1 64.1425 5.3681 -5.972 43.396 2 -2.5333 7.5098 -0.429 1.275 3 -2.8660 -5.5862 -1.631 3.095 4 1.4960 -0.1203 0.225 1.158
5 -0.6244 1.2777 0 0
(
0)
Re 4
R R
R
−
= Φ
υλ
pomiar składał się z trzech faz (rys. 4): pomiaru siły normalnej, pomiaru siły stycznej oraz tary. Podczas pomiaru siły normalnej przy zerowym kącie pochylenia oś wagi pokrywała się z kierunkiem przepływu czynnika tunelu. W takiej konfiguracji dokonywany był pomiar tarujący wagę (pomiar wypadkowych sił ciężkości i wyporu oraz pochodzących od nich momentów) dla różnych kątów pochylenia. Do pomiaru siły stycznej waga była obracana o 90o względem głównej osi robota.
System pomiarowy dokonywał pomiarów z częstotliwo- ścią 50 Hz. Sygnał dodatkowo był podawany obróbce filtrem eliptycznym o częstotliwości 100kHz. Podczas każdego pomiaru robot wykonywał trzydzieści powtó- rzeń ruchu trzepoczącego. Wyniki pierwszych trzech były odrzucane, pozostałe uśredniano. W wyniku pomiaru uzyskiwane były spadki napięć na mostkach tensometrycznych. Aby uzyskać wynik w postaci sił, należało przemnożyć spadki napięć przez odwróconą macierz kalibracji. Po przemnożeniu i uśrednieniu wyniki były aproksymowane szeregiem Fouriera o długości 20 wyrazów. Zabieg ten ułatwia dalsze przetwarzanie danych oraz dodatkowo wygładza wyni- ki. Oczywiście w uzyskanym w ten sposób wyniku są siły zarówno hydrodynamiczne, jak również siła ciężko- ści i wyporu oraz bezwładności (te właściwie są pomi- jalnie małe).
Rys. 4. Sposób pomiaru dwóch składowych sił
Aby wyodrębnić siły hydrodynamiczne, dokonywana była tara, czyli robot był ważony w bezruchu. Wobec tego siłę hydrodynamiczną można wyrazić w postaci różnicy:
I tara zmierz
hyd
F F F
F = − −
, (4)gdzie FI to siły bezwładności. Siły bezwładności zostały wyliczone na podstawie znajomości geometrii modelu oraz sposobu jego ruchu, ponieważ pomiar tych sił byłby kłopotliwy. Pomiar musiałby być wykonany poza
tunelem, więc z powodu braku siły wyporu przetwornik byłby bardziej obciążony (mógłby nawet zostać przecią- żony), co w dalszej kolejności mogłoby prowadzić do innych wyników (w wadze wykorzystane zostały tenso- metry piezoelektryczne, które mają nieliniową charakte- rystykę poza punktem pracy). Wobec powyższego zdecydowano się na prostszy sposób. Składowa normal- na siły bezwładności wyliczona została z zależności:
( ) ( ( ) t ) m r ( ) t ( ( ) t )
m I
t
F I
s scs
In
ε γ ω γ
γ γ γ
γ
γ
⋅ ⋅ cos + ⋅ ⋅
2⋅ sin
=
, (5)
natomiast składowa styczna:
(6)
Dodatkowo wyniki eksperymentu zostały porównane z wynikami obliczeń quasi - stacjonarnych, gdyż chciano się w ten sposób dowiedzieć, w jakim stopniu uproszczony model aerodynamiczny odbiega od rzeczy- wistości. Ponadto obliczone wyniki stanowią punkt odniesienia dla wyników pomiaru. Obliczenia bazowały na eksperymentalnych charakterystykach profilu (pła- skiej płytki). Do wyznaczenia charakterystyki skrzydła posłużono się prostą metodą siatki wirowej, w której pojedyncza sekcja reprezentowana jest przez pojedyn- czy element wiru. Do uwzględnienia generowanej siły podczas zmiany kąta natarcia skorzystano z teorii Theodorsena. Opis zbliżonej metody obliczeniowej można znaleźć w [8] i [9]. Należy zaznaczyć, że przyjęty model nie odzwierciedla niestacjonarnej natury opływu (nie odwzorowuje tworów wirowych śladu aerodyna- micznego, opływ z separacją odzwierciedlony poprzez zadaną funkcję opisującą współczynniki siły nośnej oporu sekcji).
5. WYNIKI EKSPERYMENTU
Na rys. 5 zestawiono przetworzone wyniki pomiaru dla zawisu (prędkość postępowa równa 0m/s). Są to prze- biegi sił normalnej do średniej płaszczyzny trzepotania (N) oraz stycznej (T). Zamieszczono również wyniki obliczeń quasi-stacjonarnych. Analizując przebieg siły
( ) ( ( ) ) ( )
( )
( ) ( ( ) )
( )
t( ( )
t)
m r( )
t( ( )
t)
r m
t t
m I
t I t
m I
t I F
sc s sc
s
s s
It
γ ω
φ ω
φ ε γ
ε φ
γ γ φ
φ
γ γ γ φ
φ φ
cos sin
cos sin
cos
⋅
⋅
⋅ +
⋅
⋅
⋅ +
+
⋅
⋅ ⋅ +
⋅ ⋅
=
2 2
normalnej, widać dwa maksima, których położenia odpowiadają maksymalnym prędkościom obrotowym skrzydła. Ruch skrzydła do przodu (I półcykl) odbywał się z mniejszą prędkością. Wartość piku dla tej fazy jest znacznie (niemalże dwukrotnie dla I połowy 0,232N, dla II 0,41N) mniejsza. Wartość minimalną (-0,013N) siła normalna uzyskuje w 1,5s okresu. Między pierwszym, a drugim maksimum przez dłuższy czas (około 1s) utrzymuje się niemalże zerowa wartość siły. Między drugim a pierwszym wartość jest wyraźnie dodatnia, lecz również przez około 1s utrzymuje się na stałym poziomie. Przebieg siły stycznej jest podobny w analo- giczny sposób do siły normalnej. Ekstrema są w tych samych miejscach. Ich wartości również znacznie różnią się, drugie jest znacznie większe. Porównując wyniki obliczeń z tymi zmierzonymi można zauważyć, że prze- bieg siły normalnej jest dość podobny. Jest jednak kilka znaczących różnic. Uzyskana na drodze obliczeń quasi - stacjonarnych wartość maksymalna składowej normalnej siły jest nieco większa, a położenie maksimum dla I półcyklu jest przesunięte. Ponadto przebieg wyraźnie przechodzi na stronę ujemną. Z porównania przebiegu sił stycznych wynika, że maksimum tych obliczonych (II półcykl) jest dużo mniejsza.
Rysunki 6 i 7 zawierają zestawienia wyników pomiarów z prób z różnymi prędkościami postępowymi, ale dla kątów pochylenia, przy których składowa siły równo- legła do wektora prędkości wody była najbliżej zera (punkt równowagi). Na wykresy naniesiono również orientacyjne przebiegi zmian kątowego położenia azy- mutalnego skrzydła(ϕ), kata nastawienia(γ) oraz pręd-
kości kątowej
∂
= ∂ t
ω φ
(wartości nie odpowiadajązamieszczonej skali na osi wartości). Wraz ze wzrostem prędkości postępowej rosną maksymalne wartości w I fazie ruchu (ruch skrzydła do przodu), maleją natomiast w fazie drugiej (ruch skrzydła do tyłu).
Znaczące różnica jest także pomiędzy 1-2,5s. Przebiegi coraz wyraźniej przechodzą poniżej osi zero. Analizując przebiegi zmierzone dla posuwów względnych 0,48 i 0,54 zauważalny jest jedynie lekki spadek wartości maksymalnych.
Rys. 5. Zmierzone i wyliczone składowe siły hydrodynamicznej dla zawisu
Rys. 6. Składowa normalna siły aerodynamicznej dla różnych prędkości postępowych
Rys. 7. Składowa styczna siły aerodynamicznej dla różnych prędkości postępowych
W przypadku sił stycznych najbardziej odróżniającym się wynikiem jest ten zmierzony dla zawisu. On cechuje się największą wartością maksymalną i asymetrią prze- biegu, więc także wypadkowa siła jest największa.
Zmiany w wartościach ekstremów są takie jak dla siły normalnej. Również różnica pomiędzy wynikami dla największych prędkości jest taka jak w przypadku pierwszej składowej siły.
W celu łatwiejszej analizy uzyskany układ sił został przetransformowany do układu współrzędnych związa-
nych z przepływem. Tak więc siła nośna zostanie wyra- żona:
) sin(
)
cos( PA T PA
N
V = ⋅ − ⋅
, (7)a siła napędowa (równoległa do kierunku przepływu):
) sin(
)
cos( PA N PA
T
H = ⋅ + ⋅
. (8)Analizując uzyskane średnie siły nośne i napędowe w dziedzinie kąta pochylenia (rys. 8), obiektu można zauważyć, że w przebadanym zakresie ich zmiany mogą być aproksymowane funkcjami liniowymi, co jest bardzo obiecujące z punktu widzenia sterowania. Siła normalna w funkcji kąta pochylenia przejawiała ten- dencję rosnącą dla dwóch prędkości postępowych J=0,1 oraz J=0,48, w pozostałych przypadkach wraz ze wzrostem kąta pochylenia (pochylenie do przodu) średnia wartość siły malała. Siła styczna do posuwu 0,31 jest malejącaw dziedzinie kąta pochylenia w prze- dziale 0,36-0,48 wykazywała charakter rosnący, nato- miast dla największej prędkości ponownie była maleją- ca. Mimo to w każdym przypadku siła napędowa rośnie wraz z pochyleniem, oczywiście wzrost ten spowodowa- ny jest kosztem spadku siły nośnej. Średnie siły nośne uzyskane na drodze obliczeń quasi-stacjonarnychnie wykazują dobrej zgodności z eksperymentem, są dużo niższe (o około 30%). Dla każdej prędkości zwiększenie pochylenia powoduje zmniejszenie wartości siły. Różni- ca pomiędzy siłami napędowymi zmierzonymi i obliczo- nymi jest mniejsza niż w przypadku siły nośnej, lecz również są mniejsze. Kąt pochylenia, jaki jest potrzeb- ny do zapewnia zerowej siły poziomej (rys. 9), czyli równowagi, zmienia się także liniowo w funkcji prędko- ści w całym zbadanym zakresie. Nie wykazuje tenden- cji do załamania. Zaobserwowana zmiana wspomniane- go kąta jest mniejsza niż wynikałoby to z modelu quasi-stacjonarnego. Na rys. 10. przedstawiona jest zmiana siły nośnej dla kąta pochylenia odpowiadające- mu równowadze. Wszystkie pomiary przeprowadzone w ruchu postępowym dają lepsze efekty niż w zawisie.
Szczególnie duży skok jest dla najmniejszej prędkości.
W zakresie posuwu względnego 0,15-0,41 siła utrzymy- wała się na stałym poziomie. W dwóch ostatnich pomiarach nastąpił ponowny skok jej wartości.
Rys. 8. Średnie wartości składowych siły aerodynamicznej dla prędkości 0,052m/s w funkcji kąta pochylenia
Rys. 9. Kąt pochylenia dla warunku równowagi
Rys. 10. Siła nośna w warunkach ustalonych
Aby odnieść uzyskane wyniki do rzeczywistości należy w pierwszej kolejności wyniki sprowadzić do postaci bez- wymiarowych współczynników:
(9)
Wyznaczony według powyższego wzoru średni współ- czynnik siły nośnej dla maksymalnej mierzonej prędko- ści wyniesie CL=2,04. Następnie poprzez przekształce- nie tego samego wzoru i podstawienie parametrów entomoptera uzyskana zostanie wymagana częstotli-
( )
∫ ⋅ Φ ⋅ ⋅ ⋅
=
RR L
dr r f r
c C V
0
2
2)
( ρ
wość trzepotania potrzebna do uzyskania pożądanej siły nośnej (równej sile ciężkości) f=11,5Hz. Teraz można wyznaczyć prędkość postępową poprzez przekształcenie wyrażenia (3.). Największa przebadana prędkość odpo- wiada 3,17m/s lotu obiektu rzeczywistego.
6.
WNIOSKINie udało się wprost zaobserwować przekroczenia mak- symalnej prędkości przy pomocy eksperymentu. Na podstawie zmian w przebiegach sił normalnych i stycz- nych w funkcji czasu (rysunki 6. i 7.) można wniosko- wać, że jest ona niewiele wyższa od tej zadanej w eksperymencie. Zmiana siły napędowej uzyskiwana jest poprzez pochylenie entomoptera. Zależność wyma- ganego kąta pochylenia dla uzyskania równowagi w funkcji prędkości lotu w rozpatrywanym zakresie jest liniowa. Siła nośna dla każdej prędkości lotu poziomego była większa od tej dla zawisu. Wykorzystany w ekspe- rymencie sposób ruchu nie jest optymalnym. Uzyskane maksymalne wartości sił stycznej i normalnej do wagi są tego samego rzędu. Oznacza to, że siła oporu, a więc także potrzebna moc jest bardzo duża. Prawdopodobnie
udałoby się poprawić sposób ruchu tak, aby napęd był efektywniejszy poprzez modyfikacje przebiegu kąta nastawienia. Generalnie należy zminimalizować siłę oporu podczas ruchu skrzydła do przodu (I półcykl).
Wypadkowa siła napędowa zostaniew ten sposób zwięk- szona.
Porównując wyniki obliczeń quasi-stacjonarnych ekspe- rymentalnymi można wyciągnąć dwa wnioski. Użyty model nie nadaje się do precyzyjnego wyznaczania osiągów entomoptera z powodu dużych rozbieżności w średnich wartościach siły nośnej oraz przebiegu siły stycznej. Zaniżenie siły stycznej (albo oporu) spowoduje zaniżenie także zapotrzebowania na moc. Z drugiej jednak strony przebieg siły normalnej był dość zbliżony, zwłaszcza w okolicy maksimów (w chwilach gdzie pręd- kość obrotowa była duża). Model wobec tego można wykorzystać do wyznaczania obciążeń w procesie kon- struowania. Być może model będzie przydatny również do modelowania dynamiki ruchu.
Opisanie powyżej badania stanowią część projektu rozwojowego finansowanego przez Narodowe Centrum Badań i Rozwoju
Literatura
1.
Czekałowsk P., Sibilski K. : Water tunnel experimental investigation on the aerodynamic performance of flapping wings for nano air vehicles. In: 29th Applied Aerodynamics Conference AIAA HONOLULU 2011.2.
RHRC, Research water tunnel specification. El Segundo California 2009.3.
Unsteady Aerodynamics for Micro Air Vehicles AC/323(AVT-149)TP/332 2010.4.
Ol, M.V.: Unsteady low Reynolds number aerodynamics for micro air vehicles (MAVs). AFRL-RB-WP-TR-2010- 3013.5.
Jong-seob Han and Jo Won Chang Flow: Visualization and force measurement of an insect based flapping wing.In: 48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition 4 - 7 January 2010, Orlando, Florida.
6.
Raport do projektu: System monitorowania i detekcji aktywności istot żywych w pomieszczeniach zakrytych i budynkach oraz ochrony i monitorowania obiektów z latającym nanorobotem w układzie entomoptera., Praca nieopublikowana.7.
Shyy W., Lian Y., Tang J., Vheru D., Liu H.: Aerodynamics of low Reynolds number flyers. Cambridge University Press: New York, 2008.8.
Lasek M., Pietrucha J., Sibilski K., Złocka M.: Modelowanie dynamiki lotu ornitoptera. Raport merytoryczny projektu nr 9 T12C 004 18, 2003.9 Leishman J.G.: Principles of helicopter aerodynamics. Cambridge University Press: New York, 2000
