SIMULATION ZUR OPTIMIERUNG DER FLUGLEISTUNGDURCH ADAPTIVE WÖLBKLAPPEN AM BEISPIEL EINES

UNBEMANNTEN KLEINFLUGZEUGES

A. Scholtz, A. KuzolapInstitut für Luft- und Raumfahrtsysteme - Technische Universität Braunschweig

Hermann-Blenk-Str. 23, 38108 Braunschweig, Deutschland

Zusammenfassung

Die wissenschaftliche und kommerzielle Nutzung von unbemannten Flugsystemen erschließt mit ihrer größer werden-den Verbreitung immer neue Anwendungsgebiete. Mit den so permanent steigenden Anforderungen an Reichweite,Flugzeit und Nutzlastmasse wird die Entwicklung neuer leistungsfähigerer Flugzeugmuster notwendig. Das Institutfür Luft- und Raumfahrtsysteme der Technischen Universität Braunschweig hat in der Vergangenheit bereits un-bemannte Flugsysteme erfolgreich in der Forschung eingesetzt. Auf dieser Erfahrung aufbauend wird ein neuesunbemanntes Mehrzweckflugzeug für die CAROLO-Familie entwickelt. Der hierbei verfolgte Ansatz einer modularenKonfiguration beinhaltet neben der Struktur auch eine rekonfigurierbare aerodynamische Auslegung, mit der eineOptimierung der Gleitzahl für verschiedene Massen und Fluggeschwindigkeiten untersucht wird. Hierzu wurde eineSimulationssoftware erstellt, mit deren Hilfe automatisiert verschiedene Wölbklappenstellungen analysiert und dieoptimalen Ruderausschläge identifiziert werden können. Für die Umsetzung wurde eine Untersuchung zur optimiertenModellerstellung durchgeführt und verschiedene Lösungsstrategien implementiert, um eine effizientes Managementder Rechenzeit zu ermöglichen.

1. EINLEITUNG

Die Nutzung unbemannter Flugsysteme (Unmanned Air-craft Systems - UAS) hat in den letzten Jahren, nichtzuletzt für die wissenschaftliche und kommerzielle An-wendung, stark zugenommen. Dabei lässt sich feststellen,dass die Vorteile geringer Größe und niedriger Kostenhäufig mit den Problemen geringer Nutzlastkapazitätenoder eingeschränkter Reichweite und Flugzeit einherge-hen.

Das Institut für Luft- und Raumfahrtsysteme (ILR) der TUBraunschweig forscht seit 2001 an der Entwicklung vonUAS. Dabei ist eine Flugzeugfamilie um den CAROLOP200 und T200 entstanden. Mit einer stetigen Erwei-terung der Anwendungsgebiete, wie beispielsweise derLandwirtschaftsforschung, dem Pipeline Monitoring, derMeteorologie oder der Forstwirtschaft [1], erfolgt auch dieNachfrage nach größeren und leistungsfähigeren Nutz-lasten. Aus diesen Gründen wurde die Flugzeugfamiliezuletzt um ein größeres modulares Flugzeugmuster mitmaximal 25 kg Abflugmasse erweitert [2], den CAROLOP360 (siehe Abb. 1). Dem Konzept der Modularität fol-gend, sollte dieses Flugzeug nicht nur eine größtmöglicheVariabilität bezüglich der Nutzlast sondern auch des Ein-satzbereiches bereitstellen, ohne signifikant an Flugzeitund Reichweite einzubüßen.

Bei herkömmlichen Auslegungsverfahren werden Flug-zeuge nur für einen Arbeitspunkt optimiert, den sogen-

nanten Auslegungspunkt. Generell kann die Verwendungvon Steuerflächen am Tragflügel durch gezielte Anpas-sung der Profilwölbung zu einer Reduktion des Wider-standsbeiwertes und damit zu einer Optimierung der fürdie Flugleistung relevanten Gleitzahl außerhalb des Aus-legungspunktes genutzt werden. Durch die in Tabelle 1angegebenen Anforderungen an das Beispielflugzeugergibt sich ein außergewöhnlich großes Spektrum anKombinationen von verschiedenen Abflugmassen undFluggeschwindigkeiten, sodass für den P360 im Folgen-den eine solche Untersuchung durchgeführt wird.

Abbildung 1: Das UAS CAROLO P360

Die hier vorgestellte Simulationsumgebung wurde erstellt,um bereits bei zukünftigen Flugzeugauslegungen das Po-tential zur Optimierung abschätzen zu können und um imvorliegenden Fall eine Einschränkung der Möglichkeiten

Deutscher Luft- und Raumfahrtkongress 2014

1

DocumentID: 340242

vor nachfolgenden Flugversuchen zur Validierung vorneh-men zu können. Da keine dreidimensionale Simulationfür höherwertige numerische Verfahren zur Lösung desStrömungsfeldes mit zufriedenstellender Genauigkeit undvertretbarem Rechenzeitaufwand gefunden werden konn-te, wurde diese im Weiteren durch eine Kombination vonXFOIL [3] und einer Erweiterung des Programms AthenaVortex Lattice (AVL) [4] selbst erzeugt.

2. VORUNTERSUCHUNGEN

Für den neuen Flugzeugentwurf steht eine deutlich er-höhte Nutzlastkapazität zur Verfügung (siehe Tabelle 1).Das modulare Strukturkonzept erlaubt es, dass Packvolu-men gering zu halten und ermöglicht darüber hinaus eineschnelle Montage auf dem Flugfeld mit einem minimalenEinsatz an benötigtem Werkzeug. Die Struktur wurde dar-auf ausgelegt, dass auch Landungen ohne Fahrwerk aufdem Rumpfboden möglich sind. Ein abwerfbares Fahr-werk bietet dabei ein großes Potential zur Widerstands-und Gewichtsreduktion, was sich leistungsteigernd aufdie Gleitzahl auswirken würde. Aus operationellen Grün-den findet diese Option in der vorliegenden Arbeit jedochkeine weitere Beachtung.

Tabelle 1: Technische Spezifikationen CAROLO P360

Parameter Wert

Zulässige Abflugmasse 15 kg bis 25 kgMaximale Nutzlast 3 kg in NutzlastbuchtBetriebstemperatur −20 ◦C bis 45 ◦CTransportgröße längstes Bauteil < 1,6mHöchstflugdauer elektrischer Antrieb 45min

Verbrennungsmotor > 2 hReisegeschwindigkeit 20m/s bis 30m/sWindgeschwindigkeit < 15m/sMontagezeit < 15min

2.1 Aerodynamische Auslegung

In die aerodynamische Auslegung des CAROLO P360fließen die guten Erfahrungen in Bezug auf die Steuerei-genschaften im Flug aus den früheren Modellen ein. Diesführt zu einem leicht überelliptischen Flügelentwurf, derdie Vorteile des geringen induzierten Widerstands mitden Auftriebsreserven im Bereich der Flügelspitze vereint[5]. Daraus ergibt sich eine gute Kontrolle des Flugzeugsim Grenzbereich zum Strömungsabriss, sowie bei Startund Landung.

Eine wichtige Anforderung besteht darin, eine meteoro-logische Sensoreinheit zu transportieren, die auf eineFluggeschwindigkeit von V = 22m/s kalibriert ist. Aus

diesem Grund wird die genannte Fluggeschwindigkeitund eine Abflugmasse von 22 kg als Auslegungspunktverwendet [5]. Um eine Optimierung der Aerodynamikauf den Einsatzfall zu ermöglichen, wird der Tragflügelmit spannweitig sechs unabhängigen Steuerflächen aus-gestattet, die eine individuelle Wölbungsverteilung ermög-lichen. Dazu wird ein Profilstrak mit einem modifiziertenProfil auf Basis der HQ/W-Familie verwendet, die auf denEinsatz mit Wölbklappen optimiert sind.

2.2 XFLR5 Simulationsergebnisse

In einem ersten Versuch [5] wurde auf Basis einer einfa-chen XFLR5-Simulation [6] ein erster Einblick in die Wirk-samkeit der Verwölbung des Tragflügels erhalten. Es istzu beachten, dass in dieser Simulation kein Rumpf- undFahrwerkseinfluss berücksichtigt wurde und die Abflug-masse aus Gründen der Vergleichbarkeit mit den in Ka-pitel 2.3 folgenden Flugversuchen auf 18,4 kg festgelegtwurde. Die Gleitzahl E wird dabei aus den verfügbarenaerodynamischen Beiwerten für Auftrieb CA und Wider-stand CW , respektive aus den im Flugversuch ermitteltenhorizontalen sowie vertikalen Geschwindigkeitskompo-nenten u und w wie folgt berechnet

E =CACW

=u

w.(1)

Abbildung 2: XFLR5 Simulationsergebnisse ohne Rumpffür κi = −5◦, 0◦ & 5◦ aus [5]

Abbildung 2 zeigt den resultierenden Verlauf der Gleit-zahl über der Fluggeschwindigkeit für die symme-trisch und betraglich gleich ausschlagenden Wölbklap-pen κi mit i = 1 bis 3. Es wurden die Klappenwinkelκi = −5◦, 0◦ & 5◦ untersucht. Die neutrale Klappenpo-sition ist dabei die klar dominierende Konfiguration. Diestrich-punktierte Linie der positiven Verwölbung zeigt im

Deutscher Luft- und Raumfahrtkongress 2014

2

Bereich von unter V = 18m/s einen geringen Vorteil, wo-gegen die negative Wölbung in dem gewählten Geschwin-digkeitsbereich keinen positiven Effekt aufzeigt. Der feh-lende Rumpfeinfluss äußert sich in allen Variationen ineiner zu hohen Gleitzahl. Da bei der Widerstandsbestim-mung das Superpositionsprinzip Gültigkeit hat, wird durcheine Berücksichtigung die Gleitzahl zu niedrigeren Be-trägen verschoben, ohne einen Einfluss auf die Positionim Geschwindigkeitsfeld zu haben. Die Interferenzeffektemit dem Tragflügel sind dabei noch nicht beachtet.

2.3 Vorläufige Flugversuchsergebnisse

Um eine Validierung der XFLR5-Simulationsergebnissezu ermöglichen, wurden Flugversuche zur Polarenmes-sung durchgeführt. Wie in [5] beschrieben, wurde auseiner Staudruck- und Statikdruckmessung des MINC-Autopilotensystems zunächst eine Geschwindigkeitspola-re bestimmt. Diese ließ sich anschließend direkt in eineGleitzahlpolare umrechnen, wie sie für diese Fälle in Ab-bildung 3 dargestellt sind. Die Daten wurden nach demHöhenstufenverfahren erhoben, für das ein stationärerGleitflug durchgeführt werden muss und die zu jederFluggeschwindigkeit gehörende Sinkrate ermittelt wird.

Abbildung 3: Mit dem Höhenstufenverfahren ermittelteFlugversuchsergebnisse für κi = −4◦, 0◦ & 4◦ aus [5]

Die Abflugmasse des Flugzeugs lag bei den oben ange-gebenen 18,4 kg und die verwendeten Wölbklappenaus-schläge entsprachen hier κi = −4◦, 0◦ & 4◦. Die um einGrad geringeren Klappenausschläge kamen durch denanliegenden Ruderdruck im Flug zu Stande, haben aberauf den qualitativen Verlauf der Polaren keinen großenEinfluss.

Im Vergleich zu den Simulationsergebnissen fällt die Ver-teilung der Gleitzahlmaxima der einzelnen Konfiguratio-nen über die Geschwindigkeit auf. So steigert die po-

sitive Wölbung die Gleitzahl im unteren Geschwindig-keitsbereich und die negative Wölbung zeigt Vorteile beierhöhter Fluggeschwindigkeit. Auch die Gewinne gegen-über der neutralen Klappenstellung fallen sehr deutlichaus, sodass die Wirksamkeit der Profilverwölbung signi-fikant erscheint. Die Höhe der Gleitzahl ist im Vergleichzur Simulation, wie vermutet, geringer. Einschränkendmuss darauf hingewiesen werden, dass die verwendeteMesssensorik keine Korrektur der äußeren Umgebungs-bedingungen zulässt und die Versuche manuell erflogenwerden müssen. Ebenso wären zur statistischen Absiche-rung mehrere Flüge mit der gleichen Klappenstellung nö-tig. Die Erkenntnisse aus diesen Versuchen rechtfertigendennoch eine umfassendere Untersuchung mit anderenVerfahren, wie sie im Folgenden vorgestellt werden.

3. DIE BEIWERTSIMULATION

Die Berechnung der Gleitzahl erfordert die Kenntnis derFlugzeugumströmung, die heutzutage mit diversen Me-thoden ermittelt werden kann. Für die Wahl der Berech-nungsmethode ist die Genauigkeit der Modellierung phy-sikalischer Effekte in der Strömungsmechanik, als auchdie Rechenzeit ausschlaggebend. Um eine durch einProgramm gesteuerte widerstandsarme Wölbklappen-stellung zu identifizieren, sind mehrere Rechnungen not-wendig. Da sich zusätzlich ein ganzes Feld von Flug-zuständen für unterschiedliche Fluggeschwindigkeitenund Gesamtabflugmassen ergibt, wird für die Umsetzungdes automatisierten Programms das Wirbelleiterverfah-ren, im folgenden als Vortex Lattice Method (VLM) be-zeichnet, verwendet. Die VLM ermöglicht eine schnelleund zuverlässige Berechnung von Beiwerten, die nähe-rungsweise gute Übereinstimmungen für Flugzuständemit kleinen Anstell-, und Schiebewinkeln, sowie geringenFluggeschwindigkeiten liefert. Im Gegensatz zu höher-wertigen Verfahren, die die Navier-Stokes Gleichungenmit verschiedenen Diskretisierungsmethoden und Turbu-lenzmodellen numerisch lösen, ist die VLM für eine unbe-aufsichtigte Anwendung gut automatisierbar und benötigtdeutlich weniger Rechenzeit.

Die Berechnung mit der VLM basiert auf Lösungen derPotentialtheorie [7], mit der der Auftrieb und der induzierteWiderstand ermittelt werden können. Es ist zu beachten,dass die Potentialtheorie eine reibungs- und drehungs-freie Strömung voraussetzt. Zusätzlich wird die Dicke desProfils in der VLM vernachlässigt und die Skelettlinie ver-wendet, da die Profildicke für die Auftriebsverteilung uner-heblich ist. Dabei wird die Tragflügelfläche mit viereckigenFlächenelementen diskretisiert und jedes Element miteinem Hufeisenwirbel belegt. Dieser Hufeisenwirbel be-steht aus einer beidseitig endlichen tragenden Wirbelliniemit der Zirkulation Γ, die an der l

4 -Linie des Elementesangeordnet wird. An deren Ende ist jeweils eine semi-

Deutscher Luft- und Raumfahrtkongress 2014

3

unendliche nicht tragende Wirbellinie angeschlossen, dieam Ende des Flächenelementes stromabwärts ausgerich-tet ist und die gleiche Zirkulation wie die gebundene Wir-bellinie besitzt [7]. Eine Berechnung des Druckfeldes mitden einzelnen Zirkulationen vereinfacht sich mit Hilfe derSuperposition, welche aufgrund der Eigenschaften derlinearen Laplaceschen Potentialgleichung angewendetwerden kann [8]. Mit der kinematischen Strömungsbedin-gung und dem Biot-Savartschen Gesetz wird ein linearesGleichungssystem, welches mit numerischen Verfahrengelöst werden kann, für die unbekannten Zirkulationenaufgestellt. Der Satz von Kutta-Joukowsky ermöglicht dieBerechnung eines Auftriebbeiwertes zu jedem tragendenWirbel anhand der Zirkulation [7]. Der induzierte Wider-stand ergibt sich aus der spannweitigen Auftriebsvertei-lung und Abwindgeschwindigkeit. Die induzierte Abwind-geschwindigkeit für einen Elementarflügel resultiert ausden nicht tragenden Wirbellinien der einzelnen lokalenFlächenelemente und lässt sich mit der Zirkulation unddem Biot-Savartschen Gesetz ermitteln.

3.1 Auswahl eines geeigneten Programms

Es sind zahlreiche Implementationen der VLM verfügbar,wie beispielsweise Athena Vortex Lattice (AVL) [4] desMassachusetts Institute of Technology (MIT), XFLR5 [6],Tornado, Surfaces (Flight Level Engineering) und VOR-LAX, die unterschiedliche Erweiterungen und Optionenbieten.

Das VLM Programm soll aus Gründen der Automatisie-rung einerseits durch Matlab gesteuert werden könnenund andererseits möglichst genaue Ergebnisse zum Flug-versuch liefern, sodass diese im Anschluss validiert wer-den können. Da die Anwendung von XFLR5 und Surfacesüber eine grafischen Benutzeroberfläche erfolgt und kei-ne Steuerung durch andere Programmiersprachen oderder Kommandozeile geboten wird, entfallen diese beidenProgramme aus der näheren Auswahl. Der AVL Quelltextberuht auf Arbeiten von J.E. Lamar der National Aeronau-tics and Space Administration (NASA) Codes [9] und demvon E. Lan, sowie L. Miranda entwickeltem ProgrammVORLAX [4]. Da AVL fortwährend gewartet wird und Ver-besserungen, sowie Erweiterungen erfolgen, wird dasProgramm VORLAX nicht in Betracht gezogen.

In dem Handbuch Surfaces [10] ist eine Bewertung unter-schiedlicher VLM und dreidimensionaler Panelverfahrenzu finden. Als Basis für die Bewertungen dienten Berech-nungen von 12 Derivativa für die Cessna 172 mit denunterschiedlichen Programmen im Vergleich zu den ver-öffentlichten Daten der Cessna Aircraft Company, wobeieine Vielzahl von Berechnungen aus [11] stammen. Fürdie Bewertung wurde jeweils die relative Abweichungdes ermittelten Wertes zu den Daten der Cessna AircraftCompany und auf diese bezogen ermittelt. Den Abwei-

chungen der fünf Programme wurden Ganzzahlen von 5für die geringste und 1 für die größte Abweichung zuge-ordnet. Aus [10] ergab sich ein Wert von 40 für AVL und31 für Tornado, sodass für die erstellte Anwendung AVLselektiert wird.

3.2 AVL

3.2.1 Grundlagen

Da AVL auf der Potentialtheorie beruht, kann nur derinduzierte Widerstandsanteil berechnet werden, der fürein Flugzeug nur einen Teil des Gesamtwiderstandesdarstellt. Um realistische Ergebnisse erzielen zu können,wird zusätzlich der viskose Profilwiderstand ermittelt, dabeide Anteile ungefähr 90% des Gesamtwiderstandesentsprechen. Gerade bei bei größeren Wölbklappenaus-schlägen steigt der Formwiderstand erheblich an undbeeinflusst die resultierende Gleitzahl maßgeblich. Fürdie Berücksichtigung des Profilwiderstands innerhalb AVLsind zwei Funktionen verfügbar.

Es besteht die Möglichkeit einen Nullwiderstand CW0 derForm

CW = CW0 +C2A

πΛe(2)

anzugeben, wobei Λ der Streckung und e der Oswald-Zahl entspricht. Dieser Nullwiderstand wird unverändertzum induzierten Widerstand addiert und besitzt keinerleiEinfluss auf die Auftriebsverteilung oder die Berechnungdes induzierten Widerstandanteils. Dieser Nullwiderstandwird in AVL als CDvis angezeigt.

Eine weitere Methode berücksichtigt die Angabe von dreiPunkten der Lilienthalpolare. Hierbei müssen die Wertefür den maximalen und minimalen Auftrieb, sowie mini-malen Profilwiderstand für jede definierte Flügelsektionangegeben werden. AVL approximiert mit diesen Stütz-punkten die Widerstandspolare mit Hilfe einer parabo-lischen Funktion. Für den ermittelten lokalen Auftriebs-beiwert innerhalb der Flügelsektion wird ein zugehörigerviskoser Widerstandsanteil aus der zweidimensionalenPolare mit der Näherung ermittelt und zum induziertenWiderstandsanteil addiert. Wenn der maximale Auftriebüberschritten oder der minimale Auftrieb unterschrittenwird, erfolgt anhand einer in AVL hinterlegten Funktionein starker Widerstandsanstieg. Mit der potentialtheoreti-schen Betrachtung ist ein Auftriebszusammenbruch nichtermittelbar, jedoch stellt die in AVL hinterlegte Methode,den Widerstandsbeiwert verstärkt zu erhöhen, eine geeig-nete Lösung für die Ermittlung einer widerstandsarmenKlappenstellung dar.

Deutscher Luft- und Raumfahrtkongress 2014

4

3.2.2 Modifikationen in AVL

Untersuchungen ergaben, dass die in AVL näherungswei-se ermittelten viskosen Widerstände sowohl Klappenaus-schläge als auch Effekte von geringen Reynoldszahlennicht korrekt abbilden [12]. Abbildung 4 stellt exempla-risch gewählte Profilpolaren einzelner Sektionen für unter-schiedliche Flugzustände, sowie Wölbklappenausschlä-ge κ dar. Hierfür wurden Profilpolaren mit XFOIL ermitteltund mit den in AVL hinterlegten Funktionen interpoliert.Die sich ergebende Reynoldszahl wurde mit der Anström-geschwindigkeit und der jeweiligen Flügeltiefe der Sek-tion gebildet. Für die Angabe der Mach-Zahl wurde dieSchallgeschwindigkeit in 0 m Höhe der Standardatmo-sphäre zugrunde gelegt. Der kritische Amplifikationsfak-tor wurde bei dem in XFOIL hinterlegtem Standardwertvon Ncrit = 9 belassen. Eine Recherche bestätigte, dassdieser Wert in Höhen von ungefähr 900 m mit dem größ-ten anzunehmendem Turbulenzgrad und der Korrelationnach Mack [13] realistisch ist [12].

In Abbildung 4a ist eine Profilpolare bei einer von derAuslegung abweichenden höheren Geschwindigkeit von33 m/s und einem Klappenausschlag von κ = −5,0◦ dar-gestellt. Es ist eine Laminardelle erkennbar, die mit dreiPunkten und parabolischen Funktionen nur näherungs-weise beschrieben werden kann, selbst wenn die Polar-punkte einer bestimmten Auswahl unterliegen, bei denendie Abweichungen möglichst gering ausfallen. Für dieAuswahl der Polarpunkte in der Abbildung wurde die Mat-lab interne max() und min() Funktion verwendet.

Der Profilpolarenverlauf in Abbildung 4b zeigt die Effektegeringer Reynoldszahlen an der Flügelspitze mit einemKlappenausschlag von κ = 15,0◦. Für unterschiedlicheAnstellwinkel verschiebt sich die Lage der Transition so-wohl auf der Ober- als auch der Unterseite sprunghaft,sodass starke Änderungen im Widerstand auftreten. Derminimale Widerstandsbeiwert ergibt sich für einen Auf-triebsbeiwert von CA = 1,1106, welcher einem Anstell-winkel von α = 1,5◦ entspricht. Für diesen Wert liegt dieTransition auf der Ober- und Unterseite bei ca. 70 % derProfiltiefe. Dabei ist zu beachten, dass die Transitions-lage sich nicht linear und sehr stark zum Anstellwinkeländert. Es ist ersichtlich, dass trotz der Auswahl bestimm-ter Polarpunkte der Verlauf des Graphen unzureichendwiedergegeben wird und die Deviationen gerade in demfür ein Flugzeug üblichen Auftriebsbereich zu groß aus-fallen.

Um die gerade bei geringen Reynoldszahlen entstehen-den Verläufe der Polaren in der Berechnung zu berück-sichtigen, wurde eine Funktion in AVL implementiert, diees ermöglicht gesamte Profilpolaren zu übergeben. Zu-vor wurden die Profilpolaren für die spannweitige Lageder Elementarflügel linear zu den charakteristischen dreiPunkten interpoliert. Jedoch ist eine solche Interpolati-

on für unterschiedlich verlaufende Polaren nicht möglich,da beispielsweise verschiedene Profilformen für die Flü-gelsektionen möglich sein müssen. Gerade im Bereichdes Randbogens ändern sich die lokalen Reynoldszahlenaufgrund der Flügeltiefe erheblich, sodass unterschied-liche Profilverläufe nicht vermeidbar sind. Unter der An-nahme, dass die in Abbildung 4 dargestellten Polarenzwei in Spannweitenrichtung aufeinanderfolgende Ver-läufe darstellen, ist eine Interpolation zwischen den zweiPolarpunkten bei minimalem Widerstand nicht sinnvoll.Zusätzlich ist eine unterschiedliche Anzahl angegebenerPolarpunkte für die einzelnen Flügelsektionen zu beach-ten.

0 1 2 3 4 5

·10−2

−1

−0, 5

0

0, 5

1

1, 5

Widerstandsbeiwert CW

Auf

trie

bsbe

iwer

tCA

(a) Re = 7,480 · 105, κ = −5,0◦, Mach = 0,097

0 1 2 3 4 5

·10−2

−0, 5

0

0, 5

1

1, 5

2

Widerstandsbeiwert CW

Auf

trie

bsbe

iwer

tCA

AVL ApproximationXFOIL

(b) Re = 4,290 · 105, κ = 15,0◦, Mach = 0,097

Abbildung 4: Darstellung der lokalen mit AVL approximier-ten und XFOIL berechneten Widerstandspolaren

Deutscher Luft- und Raumfahrtkongress 2014

5

Es wurde eine flexible Methode gefunden die sich inSpannweitenrichtung dreidimensional ergebende Kurveaus Polaren zu interpolieren, sodass unterschiedlicheVerläufe und Polarpunkte ordnungsgemäß wiedergege-ben werden. Es wird die maximal angegebene Anzahlvon Polarpunkten für das gesamte AVL Modell ermittelt,welche verwendet wird, um eine neue Lage von Polar-punkten in jeder Sektion mit der Akima Interpolation zuberechnen. Die Methode nach Akima [14] ermöglicht einekubische Spline Interpolation, die frei von Oszillationenist. Dabei wird der globale Einfluss von gegebenen Stütz-stellen gemieden und mit einer heuristischen Funktion dieSteigung an den Übergängen anhand von fünf Punktengeschätzt. Für Profilpolaren einzelner Elementarflügelzwischen zwei Sektionen werden die einzelnen Polar-punkte abhängig von der Lage linear interpoliert.

4. MODELLERSTELLUNG

In Abbildung 5 ist das AVL Modell des CAROLO P360dargestellt. Es wurde die Annahme getroffen, dass dieLeitwerksanbindung keinen merklichen Einfluss auf dieAuftriebsverteilung am Höhenleitwerk ausübt. Aus die-sem Grund wurde sie in der Modellierung vernachlässigt.Da in AVL die Berücksichtigung von zusätzlichen Elemen-ten, wie das starre Fahrwerk oder die für Forschungszwe-cke angebrachten Messinstrumente, nicht vorgesehen ist,werden diese im Modell nicht modelliert. Ihr Einfluss be-schränkt sich auf den parasitären Widerstand und würdedie Polare hin zu geringeren Gleitzahlen verschieben.

Abbildung 5: Darstellung des CAROLO P360 Modells

4.1 Rumpf

Bei Flügel-Rumpf Anordnungen findet eine wechselseiti-ge aerodynamische Interferenz statt, die die Auftriebsver-teilung des Tragflügels im Bereich des Rumpfes signifi-kant beeinflusst. Durch den Rumpf werden am TragflügelGeschwindigkeiten induziert, die einem Abwind entspre-chen und den Auftrieb vermindern. Eine Analyse für denAnstellwinkel α = 10,0◦ mit einem zum CAROLO P360vergleichbaren Modell, bei dem der Rumpfdurchmesser,

sowie die Flügelhoch- und Rücklage übereinstimmt, er-gab eine Änderung der Gleitzahl E von −5,27 % bezo-gen auf den reinen Flügel [12]. Weitere Vergleiche vonFlügel-Rumpf Einflüssen zu experimentell ermittelten Auf-triebsverteilungen [15] zeigten, dass der Einfluss desRumpfes innerhalb von AVL aufgrund der potentialtheo-retischen Betrachtung geringer ausfällt. Da AVL nur ro-tationssymmetrische Rümpfe berücksichtigen kann, beidenen die Dicke eines Profils um die Skelettlinie rotiertwird [4], wurde zu jeder Rumpfposition ein Radius be-stimmt, der durch Rotation der Querschnittsfläche desCAROLO P360 entspricht, um die Verdrängungswirkungdes Rumpfes bestmöglich wiederzugeben.

4.2 Tragflügel

In Abbildung 6 ist die Diskretisierung des Flügels mitFlächenelementen dargestellt, welche individuell mit Wir-bellinien belegt werden. Diese Flächenelemente werdenauch als Panels bezeichnet. In [4] und [6] sind einige Ein-schränkungen und Vorschriften zur Diskretisierung mitFlächenelementen bei der VLM gegeben, da die Formder viereckigen Flächenelemente, deren Anzahl, sowiedie Lage der Aufpunkte die Beiwerte und das Ergebnismaßgeblich beeinflussen. AVL bietet die Möglichkeit diePanelverteilung der Tragflügelfläche in Spannweiten- undTiefenrichtung mit einem Biasfaktor zu beeinflussen [4].

Abbildung 6: Diskretisierung der rechten Tragflügelflächedes CAROLO P360 mit Flächenelementen

Bei der Panelierung wurde darauf geachtet, dass eineVerdichtung der Panels an den Stellen erfolgt, an denenhohe Druckgradienten, respektive Änderungen der Zir-kulationsstärke auftreten. Diese Basis stellt eine genaueErfassung der Auftriebsverteilung sicher und beschleu-nigt die Rechenzeit, da mit einer geringeren Panelanzahlbessere Ergebnisse erzielt werden. Es muss beachtetwerden, dass eine minimale Panelgröße existiert [6]. Wirddiese unterschritten liefert die VLM fehlerhafte Resultateoder führt sogar zu einer instabilen Berechnung. DieserEffekt entsteht aufgrund der mathematischen Beschrei-bung der Wirbellinien, die Singularitäten entsprechen,

Deutscher Luft- und Raumfahrtkongress 2014

6

sodass ohne genügend großen Abstand des Aufpunk-tes zur Wirbellinie unendlich große Geschwindigkeiteninduziert werden können.

Die größten Druckgradienten am Tragflügel ergeben sichin der Nähe der Steuerflächen, an der Flügelvorderkan-te und an den Flügelenden. In Tiefenrichtung wurdeals Parameter Cspace = -1.3 angegeben, welcher einerMischung von Kosinus- und Minus-Sinusverteilung ent-spricht, sodass etwas mehr Flächenelemente an der Steu-erfläche als an der Vorderkante angeordnet werden. Dieeinzelnen Verteilungen von Flächenelementen und mög-lichen Parametern sind in [4] ausführlich erläutert.

0

20

40

60

t = 19, 2s

t = 13, 4s

Rec

henz

eitt

[s]

0 10 20 30 40 50

0

1

2

3

4

5

6

N = 33, 8

Anzahl Flächenelemente in Tiefenrichtung [−]

Abw

eich

ung

[%]

Rechenzeit cA CSpace = +1.0CA CSpace = +1.0 cA CSpace = -1.3CA CSpace = -1.3

Abbildung 7: Abweichung des lokalen und Gesamtauf-triebsbeiwertes mit unterschiedlicher Anzahl von Flächen-elementen, sowie Flächenverteilungen in Tiefenrichtung

Abbildung 7 zeigt die Abweichung bei einer unterschied-lichen Anzahl von Flächenelementen, sowie Flächen-verteilungen in Tiefenrichtung für einen Rechteckflügelmit einem Wölbklappenausschlag von κ = 15,0◦ beieinem Anstellwinkel von α = 0,0◦. Die Abweichungensind auf Auftriebsbeiwerte mit 50 Flächenelementen prospannweitigen Flügelsegment und einer Kosinusvertei-lung (CSpace = +1.0) bezogen. Die AVL Berechnung er-folgte reibungsfrei. Man kann erkennen, dass bei gleicherPanelanzahl und einer geschickten Verteilung dieser, so-wohl für den lokalen Auftriebsbeiwert cA, als auch fürden Gesamtauftriebsbeiwert CA, eine geringere Abwei-chung zur Referenz erzielt werden kann. Abbildung 7visualisiert, dass eine Panelverteilung mit CSpace = -1,3gegenüber der Kosinusverteilung einen rechnerischenVorteil von 3,8 Panels bei gleicher Genauigkeit und einer30,2-prozentigen Reduktion der Rechenzeit bietet. DieserVorteil gewinnt dadurch noch mehr an Bedeutung, da erin jedem spannweitigen Panelsegment zu tragen kommt.Durch eine höhere Anzahl von Panels resultiert ein größe-

res Gleichungssystem, wodurch ein progressiver Verlauffür die benötigte Rechenzeit festgestellt werden kann.Analoge Untersuchungen wurden für die spannweitigePanelverteilung am Flügel des CAROLO P360 durchge-führt und zeigten einen ähnlichen Rechenzeitanstieg mitsteigender Panelanzahl [12].

Da der vollständige Lösungsraum für die Abflugmassevon 18,26 kg ermittelt werden soll, musste ein Kompro-miss zwischen der Rechenzeit und der erzielbaren Ge-nauigkeit getroffen werden. Für das AVL Modell wurden25 Panels in Tiefenrichtung und 58 in Spannweitenrich-tung gewählt.

5. DAS STEUERUNGSPROGRAMM -IDEAL GLIDE RATIO ALGORITHM (IGRA)

Das Hauptprogramm besteht aus separaten Matlab Skrip-ten, die alle nötigen Prozesse beinhalten, um eine idealeWölbklappenstellung zu ermitteln. Die sequentiell aus-führbaren Skripte ermöglichen dabei eine frühzeitige Kon-trolle der Zwischenlösungen auf Konsistenz und helfenso fehlerhafte Rechnungen zu vermeiden. Prinzipiell las-sen sich die einzelnen Skripte ohne Änderungen in einProgramm zusammenführen. Der Programmablauf ist inAbbildung 8 dargestellt.

Das Programm kann über eine Parameterdatei eingestelltwerden, die physikalische Referenzgrößen, XFOIL undAVL Programmeinstellungen, Pfade der Modell- und Bi-närdateien, sowie Geometrieeinstellungen beinhaltet. Esbesteht die Möglichkeit eine Bandbreite von Geschwindig-keiten und Abflugmassen anzugegeben, für die jeweilseine ideale Wölbklappenstellung ermittelt werden soll.Zusätzlich besteht die Möglichkeit die Permutation vonallen Wölbklappenausschlägen zu berechnen, um dengesamten Lösungsraum zu erfassen.

Im ersten Schritt folgt eine Vorausberechnung des An-stellwinkelbereichs der Profilpolaren für die geringste undgrößte Reynoldszahl. Diese Funktion dient dazu, den ge-samten Bereich der Profilpolaren für alle Profiltiefen undFluggeschwindigkeiten zu erfassen. Im Anschluss erfolgteine Berechnung der Profilpolaren, wobei alle Klappen-stellungen in der angegeben Bandbreite und des zugehö-rigen Inkrementes parallel berechnet werden. Zusätzlichbesteht die Möglichkeit die Berechnung manuell anhanddes anzugebenden Geschwindigkeitsbereichs zu sepa-rieren, um die Berechnung auf mehrere Computer zuverteilen und die Ergebnisse im Anschluss zu fusionieren.Nach dem Import der in XFOIL berechneten Profilpolarenin Matlab bereitet ein Programm diese auf, sodass sieeine Eineindeutigkeit im Auftriebsbeiwert besitzen. DieseProzedur ist sehr wichtig, da einerseits die in Fortran hin-terlegte Interpolation nach Akima eine Eineindeutigkeitvoraussetzt und andererseits innerhalb AVL ein eindeu-

Deutscher Luft- und Raumfahrtkongress 2014

7

tiger viskoser Widerstandsbeiwert aus der Profilpolarenbestimmt werden muss. Ablöseblasen oder ein Wiederan-legen der Strömung im Bereich maximaler oder minimalerAuftriebsbeiwerte kann innerhalb der Berechnung nichtberücksichtigt werden, da bei einer Mehrdeutigkeit wei-tere Informationen zur Verfügung stehen müssen. Dabeiexistiert kein Zusammenhang zwischen dem reibungs-frei mit AVL ermittelten Auftriebsbeiwert, der außerhalbder Profilpolaren liegt und dem mit XFOIL ermittelten,zweidimensionalen und viskosen Auftriebsbeiwert selber.Einige Profilpolaren können durch das Programm nichtautomatisiert korrigiert werden, sodass anhand einer me-nügeführten Kommandozeile und zugehörigen Grapheneine manuelle Korrektur durchgeführt werden muss.

IGRA

Permutationoder

Simplex

ProgrammEinstel-lungen

m, v, κi

AVL ModellAVL ModellTemplate

AVL

CACW = CW,vis + CW,ind

E = CA

CW

Bestim-mung

max(E)

XfoilProfile

Aufbereitungder Pro-filpolaren

Datensatz2D Profil-polaren

Eingabe fürSimplex

Abbildung 8: Programmablauf von IGRA

In dem anschließenden Schritt erfolgt die AVL Simulationfür unterschiedliche Abflugmassen und Fluggeschwindig-keiten, wobei ein stationärer Horizontalflug angenommenwird, bei dem das Nickmoment in einer Trimmrechnungdurch Höhenruderausschläge kompensiert wird. Die Be-rechnung der Permutation aller Klappenstellungen erfolgtdabei parallel, da abhängig von der Genauigkeit des Klap-penausschlages mehrere Millionen Rechnungen resul-tieren. Dieses Vorgehen deckt zwar alle Lösungen abund bietet ein sicheres Ergebnis in Bezug auf eine opti-male Klappenstellung, jedoch ist dieses Verfahren sehr

rechenintensiv und kann erst durch Nutzung mehrererRechner bewältigt werden. Die Ermittlung der optimalenGleitzahl stellt ein Optimierungsproblem mit den Einga-beparametern Masse, Geschwindigkeit und Ausschlagder Wölbklappen dar. Zur Lösung dieses Problems wur-de ein Simplex Algorithmus verwendet. Mit diesem Opti-mierungsverfahren der Numerik konnte die Rechenzeitdeutlich verkürzt werden, wobei nur 3118 Rechnungenmit der Genauigkeit von 0,5◦ bei 20 Geschwindigkeitenfür die Abflugmasse von 18,26 kg benötigt wurden. DasProgramm wurde entsprechend so implementiert, dassder Simplex ebenfalls parallel gestartet wird. Dieses Ver-fahren lieferte für die gewählte Abflugmasse im vorge-gebenen Geschwindigkeitsbereich im Vergleich zum ge-samten Lösungsraum 17 von 20 passende Gleitzahlen,wobei die mit dem Simplex-Suchverfahren ermittelten,abweichenden Gleitzahlen um 0,05% geringer ausfielen.Zur Erzielung sichere Resultate ohne Abweichungen, un-tersucht das Programm mit einem gewählten Inkrementalle möglichen Klappenstellungen im Umfeld der Lösungdes Simplex. Mit dieser Methode konnten alle Ergebnisseder Permutationsrechnung bestätigt werden.

6. ERGEBNISSE

Die resultierenden Gleitzahlpolaren sind für Abflugmas-sen von 18,26 kg und 25,00 kg in Abbildung 9 dargestellt.Dabei basieren die hier gezeigten Ergebnisse der op-timalen Gleitzahl auf Gleichklappenausschlägen κi, daUntersuchungen in [12] für differenzierte symmetrischeRuderwinkel nur marginale Verbesserungen im Promille-bereich ergaben.

20 25 30 3510

15

20

25

30

Geschwindigkeit V [m/s]

Gle

itzah

lE=

CA

CW

[−]

m = 18,26 kg κi = 0,0◦

m = 18,26 kg κi für bestes Gleiten

m = 25,00 kg κi = 0,0◦

m = 25,00 kg κi für bestes Gleiten

Abbildung 9: Optimierte Gleitzahlverläufe für zwei ver-schiedene Abflugmassen

Deutscher Luft- und Raumfahrtkongress 2014

8

Es ist ersichtlich, dass für beide Massen die maximaleGleitzahl durch Verwendung der Wölbklappen in einemgroßen Bereich der Fluggeschwindigkeiten erhöht wer-den kann. Die größten Gewinne sind hierbei im Schnell-flug zu verzeichnen, da die Grundwölbung des gewähltenProfils auf die relativ niedrige Auslegungsgeschwindigkeitangepasst ausgewählt wurde. Dies wird dadurch belegt,dass für die hohe Masse im Bereich zwischen 19− 23m/sdie neutrale Klappenstellung die beste Gleitzahl ermög-licht.

Der Einfluss der Masse auf den Verlauf der Gleitzahlpo-laren äußert sich darin, dass bei geringeren Geschwin-digkeiten mit steigender Masse die Gleitzahl von 22 auf15 abnimmt, wobei sich die Gleitzahlpunkte bei neutra-ler Klappenstellung sogar halbieren. Für höhere Flug-geschwindigkeiten kehrt sich der Verlauf ab einem be-stimmten Punkt um, der in diesem Vergleich bei 23m/sliegt. Bei maximaler simulierter Fluggeschwindigkeit steigtdie optimierte Gleitzahl damit von 24 auf 28. Eine ähnli-che Differenz ergibt sich für die Grundkonfiguration mitκi = 0,0◦.

Die Begründung für diesen Effekt ist auf die Umströmungvon Körpern und dem daraus resultierenden Auftriebsbei-wert

CA =A

ρ2V

2S(3)

zurückzuführen. Bei stationärem Horizontalflug muss derAuftrieb A gleich der Gewichtskraft G sein. Um diesenFlugzustand bei einer Geschwindigkeitserhöhung ∆V ,sowie konstanter Flächenbelastung A/S und Flughöhe,respektive Dichte ρ, aufrecht zu erhalten, muss der Ge-samtauftriebsbeiwert CA entsprechend reduziert werden.Dieser kann durch eine Änderung der Profilwölbung oderdes Anstellwinkels variiert werden, wobei ohne den Ein-satz der Wölbklappenausschläge nur die Anstellwinkel-änderung für eine konstante Geschwindigkeit verwendetwerden kann. Durch die Änderung des Anstellwinkelswird der Arbeitspunkt in der Flugzeugpolaren verscho-ben und kann durch steigende Widerstandsanteile, wiebeispielsweise des Trimmwiderstands, die Gleitzahl ent-sprechend negativ beeinflussen. Für größere Massenist eine erhöhte Geschwindigkeit von Vorteil, da in die-sem Bereich geringere Anstellwinkel erforderlich sind undder Arbeitspunkt dichter am Gleitzahlmaximum zu liegenkommt. Die Nutzung von Wölbklappenausschlägen ver-schiebt die gesamte Flugzeugpolare und nimmt auf dieseWeise Einfluss auf die Gleitzahl.

In Abbildung 10 ist auf der Abszisse die Fluggeschwin-digkeit und auf der Ordinate die Masse aufgetragen. Die-ses Feld ist mit dem Gleitzahlzuwachs von optimiertenGleichklappenausschlägen κi gegenüber der neutralenKlappenstellung in Prozent überlagert. Man kann erken-nen, dass im Langsamflugbereich bis zur Auslegungsge-schwindigkeit von V = 22m/s nur geringe Gleitzahlstei-

gerungen möglich sind. Dies entspricht den schon zuvorformulierten Erwartungen. Mit steigender Fluggeschwin-digkeit konnten mit dem Optimierungsalgorithmus IGRAanalog zu Abbildung 9 für alle Massen Gleitzahlzuwäch-se errechnet werden. Der maximale Zuwachs von 13%konnte für die geringste Masse und höchste Geschwin-digkeit bestimmt werden. Betrachtet man die in Tabelle 1definierten Reisegeschwindigkeiten, so lässt sich in derAnalyse für die höchste zulässige Geschwindigkeit einGleitzahlzuwachs zwischen 4− 8% finden.

16 18 20 22 24 26 28 30 32 34 3615

16

17

18

19

20

21

22

23

24

25

Geschwindigkeit V [m/s]

Mas

sem

[kg]

0 2 4 6 8 10 12 14

Gleitzahlzuwachs [%]

Abbildung 10: Gleitzahlzuwachs über Referenzgleitzahlfür variable Masse und Geschwindigkeitskombinationen

7. ZUSAMMENFASSUNG

In dieser Arbeit wurden zunächst die Auslegungskriterienfür das unbemannte Flugsystem CAROLO P360 vorge-stellt, die zu einer Konstruktion mit sechs über den Tragflü-gel verteilten unabhängigen Wölbklappen führte. In einerVoruntersuchung wurde eine erste einfache Simulationmit Flugversuchen verglichen. Die in der Analyse auf-gezeigten Abweichungen wiesen auf die Notwendigkeiteiner genaueren Simulation zur Optimierung der Gleitzahldurch adaptive Wölbklappen hin.

Hierzu musste das VLM-Programm AVL modifiziert wer-den, um viskose Widerstandspolaren zu implementieren,die eine realistische Abschätzung der Gleitzahl erst er-möglichen. Diese wurden im vorliegnden Fall mit XFOILermittelt. Im nächsten Schritt wurde auf die in der Mo-dellerstellung relevanten Parameter der Panelierung ein-gegangen und eine daraus folgende Abschätzung der

Deutscher Luft- und Raumfahrtkongress 2014

9

Rechenzeit durchgeführt. Die Vielzahl an verschiedenenmöglichen Modellpermutationen führte zu der automati-sierten Steuerung durch das Matlab-basierte Steuerungs-programm IGRA. Darüber hinaus wurde mit einem Sim-plex Algorithmus eine Suchstrategie implementiert, mitder die optimalen Gleitzahlen in Abhängigkeit der Wölb-klappenausschläge schnell und zuverlässig identifiziertwerden können.

Die vorgestellten Ergebnisse der neuen Simulation deu-ten auf signifikante Verbesserungen der Gleitzahl hin undkonnten den Auslegungspunkt bestätigen. Sollten sichdiese Optimierungen in den anstehenden Flugversuchenbelegen lassen, so steht hiermit ein Verfahren zur Verfü-gung, welches auch in der bemannten Luftfahrt für Treib-stoffeinsparungen oder eine Steigerung der Reichweiteverwendet werden könnte.

LITERATUR

[1] A. Scholtz, C. Kaschwich, T. Krüger, K. Kufieta,P. Schnetter, C.-S. Wilkens, and P. Voersmann. De-velopment of a New Multi-Purpose UAS for ScientificApplication. In Proceedings of the International Con-ference on Unmanned Aerial Vehicle in Geomatics(UAV-g), Zürich, September 2011, ISSN 1682-1777Vol XXXVIII-1/C22.

[2] A. Scholtz. Design and Construction of a UAV-Prototype with Emergency Landing System. Diplom-arbeit, Technische Universität Braunschweig, Institutfür Luft- und Raumfahrtsysteme, July 2009.

[3] M. Drela and H. Youngren. XFOIL 6.9 User Primer,November 2001.

[4] M. Drela and H. Youngren. AVL 3.30 User Primer,August 2010.

[5] A. Scholtz, K. Kufieta, and Voersmann P. Optimi-zation of a New Multi-Purpose UAS for ScientificApplications Using Aerodynamic Reconfiguration. In

Proceedings of the 28th Congress of the Interna-tional Council of the Aeronautical Sciences (ICAS2012), Brisbane, September 2012, paper no 949.

[6] A. Deperrois. XFLR5 - Analysis of foils and wingsoperating at low Reynolds numbers, February 2013.

[7] J. J. Bertin and M. L. Smith. Aerodynamics for engi-neers. Third Edition. Prentice Hall, 1989.

[8] E. Truckenbrodt. Lehrbuch der angewandten Fluid-mechanik, volume 249. Springer-Verlag, 1 edition,1988.

[9] R.J. Margason and J.E. Lamar. Vortex-Lattice Fort-ran program for estimating subsonic aerodynamiccharacteristics of complex planforms. Technical Re-port NASA TN D-6142, National Aeronautics andSpace Administration and Langley Research Center,Febraur 1971.

[10] Flight Level Engineering. Surfaces Vortex LatticeModule, August 2009.

[11] T. Melin. A Vortex Lattice MATLAB Implementationfor Linear Aerodynamic Wing Applications. Master’sthesis, Royal Institute of Technology (KTH), Decem-ber 2000.

[12] A. Kuzolap. Simulation der optimalen Wölbklap-penausschläge für bestes Gleiten am Beispiel desunbemannten Kleinflugzeuges Carolo P360. Stu-dienarbeit, Technische Universität Braunschweig,Institut für Luft- und Raumfahrtsysteme, 2014. InVorbeitung.

[13] L.M. Mack. Transition and laminar instability. JetPropulsion Laboratory Publication, 78, 1977.

[14] A. Hiroshi. A New Method of Interpolation andSmooth Curve Fitting Based on Local Procedures.J. ACM, 17(4):589–602, October 1970.

[15] H. Schlichting and E. Truckenbrodt. Aerodynamikdes Flugzeugs Band II., volume 514. Springer-Verlag, 3 edition, 2001.

Deutscher Luft- und Raumfahrtkongress 2014

10