Berechnung der Schallabstrahlung umstr¨omter K ¨orper auf ...

Berechnung der Schallabstrahlung umstromter Korper

auf Basis zweidimensionaler Stromungssimulationen

Studienarbeit

Bjorn GreschnerMatrikelnummer 177392

Betreuer:

Prof. Dr.-Ing. Frank ThieleHermann-Fottinger-Institut fur Stromungsmechanik

Technische Universitat BerlinAugust 2002

2 INHALTSVERZEICHNIS

Inhaltsverzeichnis

1 Einleitung 1

1.1 Problemstellung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2 Zielsetzung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

2 Stromungsmechanische Grundlagen 3

2.1 Numerische Umsetzung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2.1.1 Finite-Volumen-Methode . . . . . . . . . . . . . . . . . . . . . . . . . 3

2.1.2 Zeitdiskretisierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.1.3 Konvektionsschema . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.1.4 Turbulenzmodellierung . . . . . . . . . . . . . . . . . . . . . . . . . . 5

3 Stromungsakustische Grundlagen 6

3.1 Allgemeine Akustik . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

3.2 Die inhomogene Wellengleichung . . . . . . . . . . . . . . . . . . . . . . . . . 7

3.2.1 Konzept der Greenschen Funktion . . . . . . . . . . . . . . . . . . . . 9

3.2.2 Quelltypen in der Akustik . . . . . . . . . . . . . . . . . . . . . . . . 10

3.3 Die erweiterte Lighthill-Gleichung . . . . . . . . . . . . . . . . . . . . . . . . 12

3.3.1 Darstellung von bewegten Flachen . . . . . . . . . . . . . . . . . . . . 12

3.3.2 Die Lighthill-Gleichung mit Berandung . . . . . . . . . . . . . . . . . 13

3.4 Die Gleichung von Ffowcs-Williams und Hawkings . . . . . . . . . . . . . . . 18

3.4.1 Der FWH-Algorithmus auf einer Korperoberflache . . . . . . . . . . . 24

3.5 Zweidimensionale Schallquellen . . . . . . . . . . . . . . . . . . . . . . . . . 26

3.5.1 Die zweidimensionale Greensche Funktion . . . . . . . . . . . . . . . 26

3.5.2 Die harmonische zweidimensionale Schallquelle . . . . . . . . . . . . . 29

3.5.3 Zweidimensionale Greensche Funktion in der Frequenzebene . . . . . 30

4 Numerisches Verfahren 32

4.1 Zweidimensionale Implementierung des FW-H . . . . . . . . . . . . . . . . . 32

4.2 Schnittstelle zwischen Stromungssimulation und FW-H . . . . . . . . . . . . 34

4.2.1 Wahl der Koordinatensysteme beim FW-H . . . . . . . . . . . . . . . 36

4.2.2 Konvertierung der Stromungsdaten fur die Akustik . . . . . . . . . . 38

4.3 Implementierung des FW-H in C2Noise . . . . . . . . . . . . . . . . . . . . . 39

4.3.1 Umwandlung der Orts- in eine Zeitableitung . . . . . . . . . . . . . . 39

4.3.2 Die Bestimmung der retardierten Zeit . . . . . . . . . . . . . . . . . . 43

5 Zweidimensionaler δ-Puls 45

5.1 Konfiguration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

5.2 Auswertung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

5.3 Erforderliche Diskretisierung in x3-Richtung . . . . . . . . . . . . . . . . . . 50

5.4 Zusammenfassung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

INHALTSVERZEICHNIS 3

6 Zylinderumstromung 566.1 Stromungssimulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

6.1.1 Konfiguration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 566.1.2 Rechengitter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 566.1.3 Randbedingungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 576.1.4 Experimentelle Vergleichsdaten . . . . . . . . . . . . . . . . . . . . . 586.1.5 Auswertung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

6.2 Akustiksimulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 616.2.1 Experimentelle Vergleichsdaten . . . . . . . . . . . . . . . . . . . . . 616.2.2 Konfiguration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 636.2.3 Auswertung der FW-H-Simulation . . . . . . . . . . . . . . . . . . . . 646.2.4 FW-H auf der festen Korperoberflache f1 ausgewertet . . . . . . . . . 656.2.5 FW-H auf den durchlassigen Kontrollflachen fi ausgewertet . . . . . . 71

7 NACA 0015 727.1 Stromungssimulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

7.1.1 Konfiguration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 727.1.2 Rechengitter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 737.1.3 Randbedingungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 747.1.4 Auswertung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

7.2 Akustiksimulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 767.2.1 CAA-Vergleichsdaten . . . . . . . . . . . . . . . . . . . . . . . . . . . 767.2.2 Konfiguration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 777.2.3 Auswertung der FW-H-Simulation . . . . . . . . . . . . . . . . . . . . 787.2.4 FW-H auf der festen Korperoberflache f1 ausgewertet . . . . . . . . . 787.2.5 FW-H auf der durchlassigen Kontrollflache f2 ausgewertet . . . . . . 84

8 Zusammenfassung 85

1.2 Zielsetzung 1

1 Einleitung

1.1 Problemstellung

Die Vorhersage des entstehenden Schallfeldes bei der Umstromung von Fahr- und Flugzeu-gen ist von zunehmender Bedeutung. Im Zuge verscharfter Larmvorschriften ist die genaueund effiziente Vorhersage der Schallentstehung von umstromten Korpern, z.B. Flugzeugenbei Start und Landung, heute zwingend erforderlich fur die erfolgreiche Entwicklung undZulassung neuer Verkehrsmittel.

In der Vergangenheit hat sich die Simulation von Korperumstromungen zu einem stan-dardmaßig eingesetzten Verfahren entwickelt. Dabei interessierten bisher hauptsachlich Wider-stands- und Auftriebsbeiwerte oder gemittelte Geschwindigkeitsprofile. Diese konnten in aus-reichender Genauigkeit bei vertretbarem Aufwand und hinreichend hohen Reynoldszahlenerfolgreich simuliert werden.

In den letzten Jahren stieg jedoch das Interesse an der Simulation von instationaren Stromungs-vorgangen, um weitere Fortschritte bei komplexen, ablosenden Stromungen zu erreichen. Dernumerische Aufwand fur instationare Simulationen ist jedoch sehr hoch. Um den Aufwand inGrenzen zu halten, wird haufig zweidimensional gerechnet. Dies verringert die Problemgroßeum etwa zwei Großenordnungen.

Liegen die instationaren Felder fur Druck und Geschwindigkeiten vor, so kann man auchden entstandenen Schall in jedem Punkt auswerten. Von Interesse ist jedoch meist nichtder Schall in dem von der Stromungssimulation abgedeckten Gebieten, sondern der Schallim Fernfeld. Eine Ausdehnung der Stromungssimulation auf weit ausgedehnte Gebiete istaus mehreren Grunden nicht moglich: Die numerische Dissipation vorhandener CFD-Codesist vergleichsweise hoch, so daß die kleinen Schwankungen des Schalls schon in geringerEntfernung dissipiert werden. Die Ausdehnung der Stromungssimulation auf große Gebie-te erfordert eine adaquate Diskretisierung in den entfernten Teilen und wird damit sehraufwandig.

Ein Ausweg ist die Nutzung des Algorithmus von J. E. Ffowcs-Williams und D. L. Hawkings(FW-H) [5]. Der FW-H ist eine Losung der erweiterten Lighthill-Differentialgleichung, diedie Schallerzeugung durch bewegte Korper und die Schallausbreitung in bewegten Medienbeschreibt. Der FW-H modelliert den Schall im Fernfeld durch Integrale von Monopolen,Dipolen und Quadrupolen. Die Losung wird dabei nur durch analytische Umformungen ausder erweiterten Lighthill-DGL hergeleitet.

Dadurch ergibt sich die Moglichkeit, zunachst eine Stromungssimulation fur eine optimaleKonfiguration durchzufuhren und dabei zusatzlich die instationaren Werte auf der Kon-trolloberflache zu speichern. Darauf aufbauend erfolgt im zweiten Schritt die akustischeAuswertung der instationaren Daten mit Hilfe des FW-H -Algorithmus.

1.2 Zielsetzung

Im Rahmen dieser Arbeit soll eine zweidimensionale Variante des FW-H-Algorithmus miteinem vorgegebenen Verfahren zur Losung der 3D-FW-H-Gleichung verknupft und anhandeinfacher analytischer Testfalle validiert werden. Danach sollen zweidimensionale Stromungs-simulationen fur reale Korper durchgefuhrt werden, auf deren Basis die Schallausbreitung

2 1.2 Zielsetzung

zu simulieren ist. Der Vergleich der Simulation mit experimentellen und CFD/FW-H-Datenanderer Autoren wird angestrebt.Der erste Schritt ist die Definition und Implementierung der Schnittstelle zwischen Stromungs-simulation und FW-H. Im nachsten Schritt soll die zweidimensionale Variante des FW-Himplementiert werden. Das ausfuhrlichen Testen des FW-H soll Unterschiede zur dreidimen-sionalen Akustik zeigen. Dabei sollen geeignete einfache analytisch losbare Testfalle simuliertwerden. Den Abschluß bilden die Untersuchungen fur zwei reale Stromungskorper.Der erste Testfall ist ein Kreiszylinder mit dem Durchmesser von D = 0.02m bei einerReynoldszahl von Re = 84200 und als zweite Konfiguration ein NACA0015-Flugel bei einemAnstellwinkel von α = 20o und Re = 1.5 · 106. Die CFD-Simulationen sollen mit Hilfe desam HFI entwickelten Programms ELAN2 von T. Rung [10] durchgefuhrt werden.Da die Akustik in der Realitat bei ungehinderter Schallausbreitung immer dreidimensionalist, soll geklart werden, ob mit einem zweidimensionalen Ansatz Pegel und Frequenzen richtigwiedergegeben werden konnen.

2.1 Numerische Umsetzung 3

2 Stromungsmechanische Grundlagen

Die numerische Simulation laminarer und turbulenter Stromungsvorgange basieren auf derLosung der Bilanzgleichungen fur Masse, Impuls und Energie. Alle diese Gleichungen sind vonder gleichen Form und gehen von der Erhaltung der betrachteten Große aus. Der allgemeineCharakter der Bilanzgleichungen ermoglicht die Berechnung der gesuchten Großen immerauf die selbe Art. Die Gewichtung der Terme, aufgrund der unterschiedlichen physikalischenBedeutung, differiert je nach Transportgroße.Die allgemeine differentielle Form der Transportgleichung fur die Große φ lautet:

∂

∂t(%φ)

︸︷︷︸(Instationar)

+∂

∂xi

(%uiφ)︸︷︷︸

(Konvektion)

=∂

∂xi

(Γφ

∂φ

∂xi

)

︸︷︷︸(Diffusion)

+ Sφ︸︷︷︸(Quelle)

(2.1)

Dabei ist Γφ der spezifische Diffusionskoeffizient.

Im Rahmen dieser Arbeit werden alle Stromungen als inkompressibel betrachtet, da diemaximalen Geschwindigkeiten in allen Fallen gering gegen die Schallgeschwindigkeit sind.

2.1 Numerische Umsetzung

Fur die Stromungssimulationen wird der am HFI entwickelte und mehrfach validierte Stro-mungsloser ELAN2 (elliptic analysis of Navier-Stokes equation 2D) verwendet. Das zweidi-mensionale Rechengebiet wird dabei unter Verwendung der Finite-Volumen-Methode diskre-tisiert.Das Verfahren arbeitet semiblockstrukturiert auf Basis von krummlinigen Koordinaten, sodaß komplizierte Strukturen genau und effektiv diskretisiert werden konnen.Die Losung der Gleichungssysteme erfolgt implizit mit Stone’s SIP-Solver . Dabei werden dieGleichungen der Stromungsgroßen und eine Druckkorrekturgleichung sequentiell gelost. DieDruckkorrekturgleichung, auf Basis des SIMPLE-Algorithmus, garantiert das die Geschwin-digkeits- und Druckfelder die Kontinuitatsgleichung zu jedem Zeitpunkt erfullen.In ELAN2 sind mehrere Turbulenzmodelle fur RANS-Rechnungen implementiert und va-lidiert. Dabei sind 1-Gleichungsmodelle (z.B. von Spalart-Allmaras), 2-Gleichungsmodelle(k − ε, k − ω) sowie explizite algebraische Reynoldsspannunsmodelle nutzbar.

2.1.1 Finite-Volumen-Methode

Bei der Methode der Finiten-Volumen wird zunachst ein numerisches Gitter uber das Rechen-gebiet gelegt, dessen Linien den Koordinatenlinien entsprechen. Dadurch wird das Gebietin Kontrollvolumina unterteilt. Als Lage der zu jeweils einem bestimmten Kontrollvolumengehorenden Rechenpunkte wird der Mittelpunkt der einzelnen Volumina festgelegt, wie inAbbildung 1 gezeigt wird.Die partiellen Differentialgleichungen werden uber die Kontrollvolumina integriert. Durchdiese Integration entstehen Bilanzgleichungen fur den Fluß, die eine konservative Diskretisie-rung gewahrleisten. Dies bedeutet, daß alles, was aus einem Kontrollvolumen hinausstromt,

4 2.1 Numerische Umsetzung

in das benachbarte Kontrollvolumen hineinfließt. Im Gegensatz zur Finite-Differenzen-Methodeist die Finite-Volumen-Methode daher immer konservativ.

Abbildung 1: Finites Volumenelement

Zur Bezeichnung der Gitterpunkte wird die Kompaß-Notation verwendet. Der zu berechnen-de Punkt wird mit P bezeichnet, die umgebenden Punkte werden mit N wie North, S wieSouth, E wie East und W wie West bezeichnet. Die entsprechenden Kleinbuchstaben werdenzur Bezeichnung der Kontrollvolumen-Grenzflachen verwendet.

Man kann sich die Vorstellung eines Volumens bei zweidimensionaler Betrachtungsweiseleicht herleiten, indem man die Dicke der Flachen ∆Ax, ∆Ay in z–Richtung mit 1 annimmt( ⇒ V = dxdy1).

2.1.2 Zeitdiskretisierung

Die vollimplizite Zeitdiskretisierung in ELAN2 ist in jedem Fall numerisch stabil. Die Wahlder Zeitschrittweite kann allein unter physikalischen Aspekten erfolgen bzw. in dieser Arbeitabhangig von den Frequenzen, die man spater auflosen mochte.

Durch die Kopplung des instationaren an die anderen Terme ist die Konvergenzgeschwindig-keit innerhalb eines Zeitschrittes ein guter Indikator fur die adaquate zeitliche Auflosung.

ELAN2 bietet die Wahl zwischen Zeitdiskretisierung erster und zweiter Ordnung. Die Dis-kretierung zweiter Ordnunh basiert dabei auf zwei zuruckliegenden Zeitschritten.

2.1.3 Konvektionsschema

Zur Behandlung der konvektiven Terme sind neben dem Upwind-Differenzen-Schema (up-wind differerencing scheme - UDS) und dem Zentral-Differenzen-Schema (central differeren-cing scheme - CDS) Konvektionsschemata hoherer Ordnung implementiert. Dies sind limi-tierte Differenzenschemata (total variation diminishing - TVD) sowie das QUICK-Schema(quadratic upwind interpolation for convective kinematics - QUICK).

2.1 Numerische Umsetzung 5

Zur Anwendung kommt hier das QUICK-Schema mit einem Abbruchfehler dritter Ordnung.Dies wird durch Interpolation der Zellgrenzwerte auf Basis von drei Punkten erreicht, wobeiimmer zwei Punkte stromaufwarts und ein Punkt stromabwarts liegt (Abbildung 2).

W P E EEWW

u>0 W, P, E

u<0P, E, EE

e

Abbildung 2: QUICK basiert auf 3Punkten

Bei der Wahl des Konvektionsschemas ist zu beachten, daß Diffussions- und Konvektionster-me in den Navier-Stokes-Gleichungen miteinander gekoppelt sind. Die Konvektionsschematabehandeln die Diffusion unterschiedlich, zur approximierten physikalisch sinnvollen Diffusionkommt noch ein weiterer Term hinzu, allgemein als numerische Diffusion bezeichnet.Die numerische Diffusion des QUICK-Schemas ist im Gegensatz zum CDS-Konvektionsschemahoch einzuschatzen, ein Nachteil fur die Vorhersage der Ausbreitung der kleinen akustischenDruckschwankungen. Bei vielen praxisrelevanten Problemen wird trotz dieser Problematikhaufig das QUICK-Schema verwendet und bei geringem numerischen Aufwand werden guteErgebnisse in der Vorhersage der stromungsrelevanten Parameter erzielt, der Hauptgrundfur die Verwendung von QUICK im Rahmen dieser Arbeit.

2.1.4 Turbulenzmodellierung

Entscheidenden Einfluß auf die Qualitat der gewonnen Ergebnisse hat die Tubulenzmodellie-rung. Die Grundlage der Berechnung reibungsbehafteter Stromungen ist die Navier-Stokes-Gleichung. Mit dem Ansatz von Reynolds und anschließender Mittelung ergeben sich dieRANS-Geichungen (Reynolds averaged Navier Stokes equation - RANS). Ziel der statisti-schen Turbulenzmodellierung ist die Entwicklung eines Modells zur Losung des Schließungs-problems, um die Reynoldsspannungen bestimmen zu konnen.In ELAN2 sind diverse Ein- und Zweigleichungsmodelle implementiert und validiert. Verwen-det wird hier ein von T.Rung und F.Thiele (1996) [12] an der TUB entwickeltes k−ω-Modell.Das LLR-k − ω-Modell (local linear realizibility - LLR) [9] ist eine modifizierte Version desWilcox k − ω-Modells.

6 3.1 Allgemeine Akustik

3 Stromungsakustische Grundlagen

3.1 Allgemeine Akustik

Die Akustik ist die Wissenschaft des Schalls einschließlich seiner Erzeugung, Ausbreitungund Auswirkung.Eigenschaften des Schalls sind:

• Schall breitet sich als Welle aus. Die Ausbreitungsgeschwindigkeit betragt zum Beispielin Luft unter Normalbedingungen 340 m/s.

• Schallwellen bewegen das Fluid bzw. Medium um einen mittleren Zustand.

• Schwankungen der Zustandsgroßen, z.B. Druck, durch den Schall sind fast immer re-lativ klein.

• Schallwellen transportieren Energie.

• Schallwellen konnen Informationen transportieren.

Die quantitative Beschreibung des Schalls erfolgt normalerweise durch Zerlegung der Großenin einen Gleichanteil und einen Schwankungsanteil. Zunachst soll nur ein ruhendes Mediumbetrachtet werden. So ergibt sich fur Druck p, die Dichte % und die Geschwindigkeit ~u dieZerlegung:

p(~x, t) = p0 + p′(~x, t)

%(~x, t) = %0 + %′(~x, t) (3.1)

~u(~x, t) = ~u′(~x, t)

Dabei entsprechen p0 und %0 dem Druck und der Dichte in einem Ausgangszustand ohneSchall. Sie sind unabhangig von Ort ~x und der Zeit t. Der Gleichanteil der Geschwindigkeit istim ruhenden Medium null. Die Schwankungsanteile sind mit einem Strich (·)′ gekennzeichnet.Fur die qualitative Beschreibung des Schalls wurde eine logarithmische Maßeinheit, derSchalldruckpegel SPL, eingefuhrt:

SPL = 20 log10

( pRMS

2 · 10−5Pa

)(3.2)

Die logarithmische Skala gibt den großen Bereich, den die Amplituden der Schalldruck-schwankungen in der Praxis erreichen, wider. Der auftretende Wert von 2 · 10−5Pa ist eininternational festgelegter Referenzdruck, der etwa der Horschwelle des Menschen bei 1 kHzentspricht. Damit ergibt sich fur die Horschwelle ein Schalldruckpegel von 0 dB und dieSchmerzschwelle liegt bei ca 120 dB. Der RMS- oder auch Effektivwert der Druckschwan-kung ist mit

pRMS =

√p′(t)2 (3.3)

definiert. Die zeitliche Mittelung einer Funktion y(t) ist durch

3.2 Die inhomogene Wellengleichung 7

y(t) =1

T

∫ T

0

y(t)dt (3.4)

gegeben.Auf Basis der Grundgleichungen der Stromungsmechanik und des Schwankungsgroßenansat-zes (3.1) kann man die Wellengleichung herleiten. Wird in die Kontinuitatsgleichung

∂%

∂t+

∂

∂xi

(%ui) = 0 (3.5)

und die Euler-Gleichung (ohne Reibung und Volumenkrafte)

∂

∂t(%ui) +

∂

∂xj

(%uiuj) +∂p

∂xi

= 0 (3.6)

der Schwankungsgroßenansatz einsetzt und alle Terme hoherer Ordnung in den gestrichenenGroßen vernachlassigt, so erhalt man nach einigen Umformungen die Wellengleichung:

∂2%′

∂t2− ∂2p′

∂x2i

= 0 (3.7)

Schließlich kann %′ mit Hilfe der linarisierten Druck-Dichte-Beziehung

p′ = %′c2 (3.8)

durch p′ ersetzt werden. c ist die Schallgeschwindigkeit in Luft. Kurzt man den Ausdruck∂2

∂x2i

mit dem Laplaceoperator ∆ ab, so erhalt man die in der Literatur ubliche Form der

Wellengleichung:

1

c2

∂2p′

∂t2−∆p′ = 0 (3.9)

Diese beschreibt die Ausbreitung kleiner Storungen (| p′ |¿ p0), wenn sich das Medium inRuhe befindet.

3.2 Die inhomogene Wellengleichung

Die homogene Wellengleichung beschreibt nur die Ausbreitung der Schallwellen in einemMedium. Mochte man auch die Schallerzeugung durch die Differentialgleichung ausdrucken,so muß man Quellterme auf der rechten Seite der Wellengleichung hinzufugen und dann dieLosung dieser bestimmen.Um auf elegante Weise die Losung bestimmen zu konnen, wird zuerst das akustische Potentialmit der Annahme, das die Akustik reibungsfrei ist, definiert :

~u′ = gradφ

p′ = −%0∂φ

∂t(3.10)

1

c2

∂2φ

∂t2−∆φ = 0 (3.11)

8 3.2 Die inhomogene Wellengleichung

Statt der Wellengleichung fur den Druck kann Gleichung (3.11) gelost werden und der Druckdann anschließend mit Gleichung (3.10) berechnet werden. Auf diesem Wege kann auchgleich die Schnelle (⇔ Geschwindigkeitsschwankung) angegeben werden, ohne erst wiederdie linearisierte Euler-Gleichung zum Umrechnen zwischen Druck und Schnelle zu bemuhen.Eine besondere weitere Eigenschaft zeichnet das akustische Potential aus. Kennt man eineLosung, so konnen durch partielles Ableiten der Losung nach den Koordinatenrichtungenoder der Zeit neue Losungen der Wellengleichung erzeugt werden.Die inhomogene Wellengleichung mit der Feldfunktion q(~x, t) auf der rechten Seite lautetdann

(1

c2

∂2

∂t2−∆

)φ = 4πq(~x, t). (3.12)

Eine kontinuierliche Quellverteilung definiert man mit Hilfe der δ-Funktion:

q(~x, t) =

∫

R3

q(~y, t)δ(~x− ~y)d3~y (3.13)

Die kontinuierliche Quellverteilung q(~x, t), die im gesamten Bereich der Differentialgleichungdefiniert ist, wird so durch eine beschrankte Anzahl von Punktquellen q(~y, t) dargestellt.Dies wird erreicht durch zwei besondere Eigenschaften der δ-Funktion. Ihr Wert ist gegebendurch

δ(x) =

{1 fur x = 00 sonst

(3.14)

und das Intgral uber δ ist eins. Mit der δ-Funktion ist es somit moglich, das Integral ubereine beliebige stetige Funktion zu bilden und als Ergebnis nur den Wert der Funktion an derStelle x0 zu erhalten. Dies wird durch Multiplikation der Funktion mit δ(x− x0) erreicht.Das Feld

φ(~x, t) =

∫

R3

q(~y, t− |~x−~y|c

)

4π | ~x− ~y | d3~y (3.15)

erfullt damit die inhomogene Wellengleichung (3.12). Die Losung kann man durch Einsetzenuberprufen und auch herleiten [4].Gleichung (3.15) gilt nur im unendlich ausgedehnten Raum und kann folgendermaßen inter-pretiert werden:

Eine Punktquelle am Ort ~y erzeugt beim Beobachter in ~x zum Zeitpunkt t ein Signal derGroße q(~y,τ)

|~x−~y| . Die Punktquelle wird dabei zur retardierten Zeit

τ = t− | ~x− ~y |c

(3.16)

ausgewertet. Wie in Abbildung 3 verdeutlicht, muß das Signal die Strecke | ~x − ~y | zuruck-legen. Dabei vergeht die Zeitspanne | ~x − ~y | /c. Diese muß vom Beobachtungszeitpunkt tabgezogen werden, um die Quellzeit zu erhalten. Die retardierte Zeit hangt damit von demBeobachtungsort, dem Quellort und der aktuellen Zeit ab.


x

y

| x−y |

Quelle

Beobachter

Abbildung 3: Erklarung der retardierten Zeit

Die Losung fur mehrere Punktquellen wird durch Summieren der Losungen fur die einzelnenQuellen gebildet.

3.2.1 Konzept der Greenschen Funktion

Abschließend soll eine spezielle Darstellung der Losung (3.15) vorgestellt werden. Das Kon-zept der Greenschen Funktion versucht die Losung der inhomogenen Wellengleichung alleinedurch ein Integral uber das Produkt aus Quellterm und einer Funktion G darzustellen.Fur eine stetige Funktion B(τ) gilt:

B(τ0) =

+∞∫

−∞

B(τ)δ(τ0 − τ)dτ (3.17)

Dabei kann τ0 frei gewahlt werden. Hier wird speziell

τ0 = t− | ~x− ~y |c

(3.18)

gewahlt und die Funktion B(τ) mit

B(τ) =q(~y, t)

| ~x− ~y | (3.19)

fur feste Werte von ~x und ~y definiert. Werden beide Definitionen in (3.17) eingesetzt, so folgtder Zusammenhang

q(~y, t− |~x−~y|c

)

| ~x− ~y | =

+∞∫

−∞

q(~y, τ)

| ~x− ~y |δ(t− | ~x− ~y |

c− τ

)dτ. (3.20)

Auf der linken Seite ergibt sich der Integrand aus Gleichung (3.15). Damit kann dieser durchdie rechte Seite ersetzt werden. Man erhalt fur die Losung der inhomogenen Wellengleichungden Ausdruck

10 3.2 Die inhomogene Wellengleichung

φ(~x, t) =

∫

R3

+∞∫

−∞

G(~x, t, ~y, τ)q(~y, τ)dτd3~y (3.21)

mit

G(~x, t, ~y, τ) =δ(t− r

c− τ)

4πrmit r =| ~x− ~y | . (3.22)

Die Losung ist nun als Integral uber den Raum und die retardierte Zeit τ gegeben. Diesespezielle Greensche Funktion (3.22) gilt naturlich nur fur den dreidimensionalen Fall im unbe-grenzten Raum. Es gibt eine anschauliche Erklarung fur die einfache Darstellung der Losungder inhomogenen Wellengleichung mit der Grenschen Funktion. Die Greensche Funktion istauch Losung der Wellengleichung. Wird als Quellterm ein kurzer Puls mit der Amplitude Aeingesetzt, so erhalt man als Losung φ die Greensche Funktion G. Das heißt der Puls kommtzur retardierten Zeit mit der Amplitude A

rbeim Beobachter an. Die Greensche Funktion ist

als Elementarwelle interpretierbar.Mussen Randbedingungen beachtet werden oder zweidimensionale Schallquellen betrachtetwerden, muß eine andere passende Greensche Funktion gesucht werden. Bei Kenntnis derneuen Greenschen Funktion ist dann auch die Losung des Problems sofort mit Gleichung(3.21) bekannt.

3.2.2 Quelltypen in der Akustik

Zur Beschreibung von beliebigen Schallquellen hat sich in der theoretischen Akustik dieNutzung von Punktquellen als vorteilhaft erwiesen. Es werden drei Basistypen zur Charak-terisierung genutzt: Monopol, Dipol und der Quadrupol. Komplizierte Schallquellen werdendurch Kombination dieser Quelltypen approximiert. Die Losungen fur die drei Quelltypenkonnen analytisch aus der Definition des akustischen Potentials hergeleitet werden. Es lassensich einfache reale Schallquellen finden, deren Losungen die gleichen Eigenschaften wie diedrei Basisquelltypen besitzen. Dies ermoglicht meist eine einfache physikalische Interpretati-on von Schallereignissen. Im folgenden wird die allgemeine inhomogene Wellengleichung furden Schalldruck betrachtet.

(1

c2

∂2

∂t2−∆

)p′ = 4πq(~x, t) (3.23)

Die punktformigen Quellterme liegen in den folgenden Betrachtungen im Koordinatenur-sprung. Die Losungen fur die Wellengleichung lauten:

Monopol:

p′(~x, t) =

{h(t− |~x|

c)

| ~x |

}(3.24)

ist das Monopolfeld als Losung fur die Punktquelle qM mit h als beliebige zeitabhangigeFunktion. h beschreibt die zeitliche Anderung der Punktquelle.


qM(~x, t) = 4πh(t)δ(~x) (3.25)

Die Losung des Monopols beschreibt die in alle Richtungen gleichmaßige Schallaus-breitung. Das Signal h breitet sich kugelformig aus.

Dipol: Der Dipol laßt sich durch eine Ortsableitung des Monopolfeldes herleiten.

p′(~x, t) =∂

∂xi

{h(t− |~x|

c)

| ~x |

}(3.26)

Dieselbe Losung ist auch durch Kombination zweier Monopole darstellbar, die gegen-phasig schwingen und sich im Koordinatenursprung in geringem Abstand zueinanderbefinden. So erklart sich auch der Name Dipol (gleich Zweipol).

Betrachtet man beispielsweise die Ableitung in x1-Richtung, so loschen sich die gegen-phasigen Signale der Monopole entlang der x2- bzw x3-Achse aus. In x1-Richtung istdas Signal am großten. Es ergibt sich die typische Richtcharakteristik der Losung inForm einer acht.

Der entsprechende Quellterm lautet:

qD(~x, t) = 4π∂

∂xi

{h(t)δ(~x)

}(3.27)

Quadrupol:

p′(~x, t) =∂2

∂xixj

{h(t− |~x|

c)

| ~x |

}(3.28)

Die Losung des Quadrupols (3.28) laßt sich durch zweifache Ortsableitung eines Mo-nopolfeldes herleiten bzw. analog zum Dipol durch Kombination aus vier einzelnenMonopolen darstellen. Der entsprechende Quellterm lautet:

qQ(~x, t) = 4π∂2

∂xixj

{h(t)δ(~x)

}(3.29)

Monopol, Dipol und Quadrupol sind Losung der Wellengleichung, entsprechen jedoch keinerphysikalisch sinnvollen Quelle.

Die Wellengleichung (3.9) wird aus der Kontinuitatsgleichung (3.5) und der Eulergleichung(3.6) hergeleitet. Die Kontinuitatsgleichung basiert auf der Erhaltung der Masse. Eine Storungder Massenerhaltung, z.B. eine zeitlich schwankende Massenquelle, wird durch einen Quell-term auf der rechten Seite der Kontinuitatsgleichung dargestellt. Auf der Basis der modifi-zierten Gleichungen erhalt man bei der Herleitung der Wellengleichung einen Massequelltermqmasse auf der rechten Seite:

qmasse(~x, t) = %0∂

∂t

{β(t)δ(~x)

}(3.30)

12 3.3 Die erweiterte Lighthill-Gleichung

Die Große %0β entspricht der zugefuhrten Masse pro Zeit und pro Volumen. Die δ-Funktionsorgt dafur, daß sich bei Integration uber ein Volumen um die Quelle herum ein endlicherMassenstrom, d.h. Masse pro Zeit, ergibt.

Die Massenquelle (3.30) hat eine ahnliche Form wie der Monopol qM (3.25) . Sie unter-scheiden sich durch die zeitliche Ableitung. Die Losung wird sich also ebenso nur durch dieAbleitung des Signals von dem des Monopolfeldes unterscheiden. Die Eigenschaften des Mo-nopols konnen somit auf die Massenquelle ubertragen werden. Der von einer Massenquelleqmasse erzeugte Schall breitet sich wie beim Monopol kugelformig aus.

Analog zur Massenquelle ist es moglich, durch Storung der Impulserhaltung einen Quelltermin der Euler-Gleichung zu erzeugen. Die Impulsquelle qimpuls erscheint ebenso als Quelltermin der Wellengleichung:

qimpuls(~x, t) = − ∂

∂xi

{fi(t)δ(~x)

}(3.31)

Der Vektor fi(t) gibt die Richtung und Starke des zugefuhrten Impulses zur Zeit t an. DieOrtsableitung der Impulsquelle ist charakteristisch fur einen Dipol qD. Die Quellterme vonImpulsquelle und Dipol haben eine identische Form. Die Schallausbreitung erfolgt gerichtet.

3.3 Die erweiterte Lighthill-Gleichung

Die bisherigen Betrachtungen zur Darstellung von Schallentstehung und Schallausbreitunggelten nur fur ruhende Medien und im freien Raum. Erzeugt eine Oberflache Schall, z.Beine Lautsprechermembran, so kann man versuchen, diesen Schall durch Definition einerbegrenzten Anzahl von Monopolen, Dipolen und Quadrupolen im Quellterm q(~y, t) der Wel-lengleichung zu modellieren, und anschließend die Losung bestimmen. Der komplizierte Teilist dabei die Bestimmung der passenden kontiuierlichen Quellverteilung. Desweiteren gibt esnur fur bestimmte Spezialfalle einfache analytische Losungen. Muß man auch Randbedin-gungen erfullen, z.B. eine schallreflektierende Wand, so ist die inhomogene Wellengleichungin der Form (3.12) nicht mehr verwendbar.

Die Idee der erweiterten Lighthill-Gleichung ist die Bestimmung einer Wellengleichung furbewegte Korper aus den nichtlinearen Stromungsgleichungen. Die Reibungsterme und Nicht-linearen Terme sollen dabei auf die rechte Seite gebracht werden, damit auf der linken Seiteder normalen Wellenausdruck erhalten bleibt.

3.3.1 Darstellung von bewegten Flachen

Das Vorhandensein von Korpern erfordert die Erfullung der Randbedingungen an der Ober-flache. Voraussetzung fur eine allgemein gultige Wellengleichung ist eine formale Darstellungvon sich bewegenden Flachen, die als Feldfunktion in die Wellengleichung einfließt.

3.3 Die erweiterte Lighthill-Gleichung 13

f(x,t) < 0

f (x, t) > 0

S

n

Abbildung 4: Beschreibung der Oberflache durch f(~x, t)

Zur Beschreibung einer bewegten Flache im Raum wird eine Hilfsfunktion f(~x, t) definiert.Gegeben sei ein sich zeitlich veranderliches Volumen V mit der Oberflache s. Wie in Abbil-dung (4) verdeutlicht, wird f(~x, t) so konstruiert, daß folgende Bedingungen erfullt sind:

f(~x, t) < 0 falls ~x innerhalb von V

f(~x, t) > 0 falls ~x außerhalb von V (3.32)

f(~x, t) = 0 auf der Oberflache S

Die Hilfsfunktion soll stetig und differenzierbar sein. Zusatzlich soll die folgende Bedingunggelten:

gradf(~x, t) 6= 0 auf der Oberflache S (3.33)

Der Gradient von f(~x, t) auf der Oberflache S zeigt damit immer in Richtung des Norma-lenvektors ~n. Der Normalenvektor kann direkt aus f(~x, t) berechnet werden:

~n(~x, t) =gradf(~x, t)

| gradf(~x, t) | (3.34)

Durch die Hilfsfunktion wird die Bewegung der Oberflache eindeutig festgelegt. Definiertman eine Normalengeschwindigkeit ~v, die Geschwindigkeit, mit der sich ein Punkt auf derOberflache S in Richtung ~n bewegt, so gilt fur f(~x, t):

∂f(~x, t)

∂t+ vi

∂f

∂xi

= 0 (3.35)

3.3.2 Die Lighthill-Gleichung mit Berandung

Existieren keine Korper im Gebiet, so kann die Losung der Wellengleichung direkt angegebenwerden. Die Quellstarkeverteilung wird mit der Greenschen Funktion multiplitziert und uberden Quellbereich integriert. Da innerhalb von Korpern die Wellengleichung nicht gilt, kannman diesen einfachen Weg bei Anwesenheit von Korpern nicht mehr nutzen.


Die erweiterte Lighthillgleichung erweitert die Gultigkeit der Wellengleichung auch auf Ober-flachen bzw. Korper, so daß die Losung wieder mit der Greenschen Funktion moglich ist.Dazu benotigt man die Heaviside-Funktion H(x):

H(x) =

{1 fur x > 00 fur x < 0

(3.36)

Mit der zuvor definierten Funktion f(~x, t) als Argument von H(x) erhalt man:

H(f(~x, t)) =

{1 fur ~x 6∈ V ; f(~x, t) > 00 fur ~x ∈ V ; f(~x, t) < 0

(3.37)

Fur die Herleitung der erweiterten Wellengleichung werden die nichtlinearen Erhaltungsglei-chungen fur Masse und Impuls benutzt. Die Kontinuitatsgleichung

∂%

∂t+

∂

∂xi

(%ui) = 0 (3.38)

wird mit H(f) multipliziert und eine frei wahlbare Konstante %0 eingefuhrt. Die zeitlicheAbleitung von %0 sei null. Es ergibt sich

{∂

∂t(%− %0) +

∂

∂xi

(%ui)

}H(f) = 0. (3.39)

H(f) kann in die Ableitungen hineingezogen werden und es gilt

∂

∂t

{(%− %0)H(f)

}+

∂

∂xi

{(%ui)H(f)

}

= (%− %0)∂

∂t

{H(f)

}+ (%ui)

∂

∂xi

{H(f)

}. (3.40)

Die Zeitableitungen der Heaviside-Funktion ergibt sich zu

∂

∂t{H(f(~x, t))} =

dH

df

∂f

∂t= δ(f)

∂f

∂t. (3.41)

Der Ausdruck auf der rechten Seite ist durch die δ- Funktion nur bei f(~x, t) = 0 , d.h. aufder Oberflache S ungleich null. Fur die Oberflache gilt die Gleichung (3.35), so daß man dieZeitableitung ∂f

∂tersetzen kann. Es folgt

∂

∂t{H(f(~x, t))} = −vi

∂f

∂xi

δ(f). (3.42)

Fur die raumliche Ableitung ergibt sich entsprechend

∂

∂xi

{H(f(~x, t))} =∂f

∂xi

δ(f). (3.43)

Setzt man die Ableitungen in Gleichung (3.40) ein, so erhalt man

∂

∂t

{(%− %0)H(f)

}+

∂

∂xi

{(%ui)H(f)

}=

{%(ui − vi) + %0vi

} ∂f

∂xi

δ(f). (3.44)


Gleichung (3.44) ist eine Erweiterung der Kontinuitatsgleichung. Innerhalb des Volumens Vist die linke und rechte Seite null. Außerhalb ergibt sich die einfache Kontinuitatsgleichugnach Gleichung (3.38). Nur auf der Oberflache S ist die rechte Seite von null verschieden. Dieerweiterte Kontinuitatsgleichung kann so im gesamten Raum angewendet werden, auch wennsich Korper im Raum befinden. Es mussen an der Korperoberflache keine Randbedingungenerfullt werden. Dafur hat man auf der rechten Seite einen zusatzlichen Term, und es mußdie Hilfsfunktion f(~x, t) so gewahlt werden, daß sich der Korper innerhalb des Volumens Vbefindet. Dabei kann die Oberflache S auch mit der Korperoberflache ubereinstimmen.Analog zur Kontinuitatsgleichung kann auch die Impulserhaltungsgleichung umgeformt wer-den. Diese enthalt im Gegensatz zur Eulergleichung (wird genutzt zur Herleitung der homo-genen Wellengleichung (3.6)) auch die Reibungseffekte. Die Impulsgleichung lautet

∂

∂t(%ui) +

∂

∂xj

(%uiuj + Pij) = 0 (3.45)

mit dem Tensor

Pij = (p− p0)δij − τij. (3.46)

Dabei ist τij der Spannungstensor, und p0 ist eine frei wahlbare Konstante, deren raumlicheAbleitung verschwindet.

δij =

{1 fur i = j0 fur i 6= j

(3.47)

ist die sogenannte Kronecker-Deltafunktion. Nach Multiplikation der Impulserhaltungsglei-chung mit H(f), Umformen der Gleichung, sodaß H(f) innerhalb der Ableitungen steht undErsetzen der Ableitungen mit Gleichungen (3.42) und (3.43) ergibt sich

∂

∂t

{(%ui)H(f)

}+

∂

∂xj

{(%uiuj + Pij)H(f)

}=

{%ui(uj − vj) + Pij

} ∂

∂xj

δ(f). (3.48)

Gleichung (3.48) gilt ebenso wie die erweitert Kontinuitatsgleichung im gesamten Raum mitKorpern, sofern f(~x, t) geeignet gewahlt wird. Aus den erweiterten Erhaltungsgleichungenkann eine Wellengleichung hergeleitet werden. Dazu wird als erstes von der Impulserhal-tungsgleichung (3.48) die Divergenz gebildet.

∂

∂t

∂

∂xi

{(%ui)H(f)

}+

∂

∂xi

∂

∂xj

{(%uiuj + Pij)H(f)

}=

∂

∂xi

{(%ui(uj − vj) + Pij)

∂

∂xj

δ(f)

}

(3.49)Die Kontinuitatsgleichung (3.44) wird partiell nach der Zeit abgeleitet:

∂2

∂t2

{(%− %0)H(f)

}+

∂

∂t

∂

∂xi

{(%ui)H(f)

}=

∂

∂t

{(%(ui − vi) + %0vi)

∂f

∂xi

δ(f)

}(3.50)

Dann zieht man Gleichung (3.49) von (3.50) ab, dabei heben sich zwei Terme auf der linkenSeite gegeneinander auf. Es ergibt sich


∂2

∂t2

{(%− %0)H(f)

}− ∂2

∂xixj

{(%uiuj + Pij)H(f)

}

=∂

∂t

{(%(ui − vi) + %0vi)

∂f

∂xi

δ(f)

}− ∂

∂xi


∂

∂xj

δ(f)

}. (3.51)

Um die Form einer Wellengleichung zu erhalten, wird auf beiden Seiten der Ausdruck

c2∆{(%− %0)H(f)

}= c2 ∂2

∂xixj

{(%− %0)δijH(f)

}(3.52)

subtrahiert und der zweite Term auf der linken Seite von (3.51) durch Addition auf die rechteSeite gebracht. c ist wiederum eine frei wahlbare Konstante, die fur sinnvolle Anwendungenspater gleich der Schallgeschwindigkeit gesetzt wird. Man erhalt schließlich die erweiterteLighthillgleichung

(∂2

∂t2− c2∆

) {(%− %0)H(f)

}=

∂2

∂xixj

{TijH(f)

}

+∂

∂t

{(%(ui − vi) + %0vi)

∂f

∂xi

δ(f)

}(3.53)

− ∂

∂xi


∂f

∂xj

δ(f)

}

mit dem Lighthillschen Spannungstensor

Tij = %uiuj + Pij − c2(%− %0)δij. (3.54)

Gleichung (3.53) stellt eine inhomogene Wellengleichung fur den Ausdruck

(%− %0) ·H(f) = %′ ·H(f) (3.55)

mit Quelltermen auf der rechten Seite dar. Damit ergibt sich die kompakte Darstellung vonGleichung (3.53) zu

(∂2

∂t2− c2∆

) {%′H(f)

}= q(~x, t). (3.56)

Sie folgt ohne Naherung aus den nichtlinearen Erhaltungsgleichungen. Alle Nichtlinearitatensind im Quellterm q(~x, t) enthalten. Die Gleichung gilt uberall, auch wenn sich Korper imRaum befinden. Diese Form der erweiterten Lighthillgleichung wurde als erstes von Ffowcs-Williams und Hawkings [5] vorgestellt. Jedoch hat Lighthill [7] schon 1951 eine vereinfachteForm hergeleitet, bei nur der erste Quellterm mit dem Lighthillschen Spannungsterm vor-kommt. Er benutzte sie zur Vorhersage von Strahltriebwerkslarm.


Die Losung der erweiterten Lighthillgleichung kann wieder mit der Greenschen Funktionbestimmt werden, d.h. sie ist wieder durch folgendes Integral darstellbar.

{%′H(f)

}(~x, t) =

1

4π

∫

R3

q(~y, t− |~x−~y|c

)

| ~x− ~y | d3~y (3.57)

18 3.4 Die Gleichung von Ffowcs-Williams und Hawkings

3.4 Die Gleichung von Ffowcs-Williams und Hawkings

Die Kenntnis der Losung in der Form (3.57) ist jedoch nicht von sehr großem Nutzen, wennman ein reales Problem hat. Die freien Konstanten in der erweiterten Lighthillgleichungmussen so gewahlt werden, das die Terme auch fur die Akustik sinnvoll sind. So wird durchDefinition der freien Konstanten %0 und p0 mit den Werten fur die Ruhedichte und Ruhe-druck, d.h. die Werte, solange kein Schall vorhanden ist, auch c als einzig sinnvolle Belegungmit der Schallgeschwindigkeit festgelegt. Damit ist die Gultigkeit der Differentialgleichungebenso eingeschrankt. Das Bezugssystem muß dann ein nicht bewegtes Koordinatensystemsein, da sonst der Wellenoperator in dieser Form nicht mehr gilt. Verdeutlicht wird dies inAbbildung 5.

Beobachter

x

x

x

y

η

1

2

f (x, t)u

Abbildung 5: Integration im korperfesten Koordinatensystem ~η

Der Flugel bewegt sich im ruhenden Fluid mit der Geschwindigkeit v nach links. Der Beob-achter befindet sich in ~x, und die Quellpunkte auf der Oberflache sind mit ~y gekennzeichnet.In seiner direkten Umgebung wird der Flugel turbulente Schwankungen im Fluid verursachen(werden durch Tij erfasst), Krafte auf das Medium ausuben (da der Flugel Auftrieb erzeugt)und das Fluid verdrangen. Diese Großen schwanken in dem ~x-Koordinatensystem zeitlich furjeden Ort , es entsteht Schall. In genugend großer Entfernung sind alle Quellterme auf derrechte Seite der DGL null. Die erweiterte inhomogene Wellengleichung wird im Unendlichenzur homogenen Wellengleichung.Damit die erweiterte Wellengleichung mit den gewahlten Parametern %0, p0 und c die Schall-ausbreitung des entstandenen Schalls korrekt beschreibt, mussen alle Großen, die in der DGLauftauchen, durch das ~x-Koordinatensystem als Bezugssystem beschrieben werden.Das großte Problem stellt dann jedoch die Hilfsfunktion f dar, die die bewegte Korper-oberflache zu jedem Zeitpunkt beschreibt. Dabei muß f(~x, t) = 0 fur alle Zeiten t im ~x-Koordinatensystem gelten und danach die Integrale der Quellterme berechnet werden. Diesfuhrt schon fur einfache Bewegungen, z.B. ein sich drehender Hubschrauberrotor, auf sehrschwer zu losende Integrale.

3.4 Die Gleichung von Ffowcs-Williams und Hawkings 19

Einen Ausweg stellten Ffwocs-Williams und Hawkings [5] 1969 vor. Sie definieren die Hilfs-funktion f in einem zweiten korperfesten Koordinatensystem ~η und fuhren auch alle In-tegrationen in diesem aus. Außerdem werden die beiden Raumintegrale, die nur auf derOberflache Werte liefern, in Oberflachenintegrale umgeformt. Fur die Herleitung muß dannjedoch vorausgesetzt werden, daß die Oberflache starr ist.Die erweiterte Lighillgleichung lautet

(∂2

∂t2− c2∆

) {%′H(f)

}=

∂2

∂xixj

{TijH(f)

}+

∂

∂t

{Qi

∂f

∂xi

δ(f)}− ∂

∂xi

{Lij

∂f

∂xj

δ(f)}

(3.58)

mit

Tij = %uiuj + Pij − c2(%− %0)δij

= %uiuj − τij + δij(p′ − c2%′)

Qi = %(ui − vi) + %0vi (3.59)

Lij = %ui(uj − vj) + Pij

Pij = (p− p0)δij − τij.

Dabei ist ui die Geschwindigkeit des Fluids und vi die Geschwindigkeit der Oberflache S, diedurch f(~x, t) = 0 beschrieben wird. Es gilt fur f auf der Oberflache

∂f

∂xi

= ni | gradf | . (3.60)

Es werden jetzt nur die beiden Terme Qi und Lij betrachtet, die auf der Oberflache Werteliefern. Der Quellterm q(~x, t) lautet dann

q(~x, t) =∂

∂t

{(%(ui − vi)ni + %0vini) | gradf | δ(f)

}

− ∂

∂xi

{(%ui(uj − vj)nj + Pijnj) | gradf | δ(f)

}. (3.61)

Fur die flachenhafte Modellquellenverteilung

q(~x, t) = Q(~x, t) | gradf | δ(f) (3.62)

und die inhomogene Wellengleichung

(1

c2

∂2

∂t2−∆

)p′ = q(~x, t) (3.63)

ergibt sich die Losung

p′(~x, t) =1

4π

∫

R3

1

rQ(~y, τ) | gradf(~y, τ) | δ(f(~y, τ))d3~y (3.64)


mit

τ(~x, ~y, t) = t− | ~x− ~y |c

und r(~x, ~y) =| ~x− ~y | . (3.65)

In Gleichung (3.64) sind die Funktionswerte Q(~y, τ) und f(~y, τ) im Integrand zur retardiertenZeit τ zu nehmen. Der Integrand ist fur ein festes ~x und t nur an den Stellen ~y mit

f(~y, t− | ~x− ~y |

c

)= [f ]ret = 0 (3.66)

von Null verschieden. Die Darstellung mit der eckigen Klammer [·]ret gibt an, daß der inihr eingeschlossene Ausdruck zur retardierten Zeit am Quellort ausgewertet werden muß.Gleichung (3.64) wird mit einem Integral uber die retardierte Zeit erweitert.

p′(~x, t) =1

4π

∫

R3

+∞∫

−∞

1

rQ | gradf | δ(f)δ(g)dτd3~y (3.67)

Jetzt ist τ eine unabhangige Integrationsvariable und entspricht damit nicht mehr der spe-ziellen retardierten Zeit aus Gleichung (3.65). Diese spezielle retardierte Zeit wird fortanmit τ ? bezeichnet. Die Variablen ~y und τ sind unabhangig, die Integrationen konnen somitvertauscht werden. Die Funktionen Q und f sind weiterhin nur abhangig von y und τ . DieFunktion g muß dann folgendermaßen definiert sein:

g(~x, t, ~y, τ) = t− | ~x− ~y |c

− τ (3.68)

Fur eine feste Beobachtungsposition ~x und Beobachtungszeit t ist die Funktion g nur von yund τ abhangig. Der Integrand ist nur an den Stellen mit f(~y, τ) = 0 und g(~y, τ) = 0 vonNull verschieden.

Fur den eindimensionalen y-Raum sind die Funktionen f = 0 und g = 0 in Abbildung 6 in der(y, τ)-Ebene dargestellt. Die Bedingungen g = 0 entspricht zwei geraden Strahlen (mit derSteigung c), die sich im Beobachtungspunkt (x, t) treffen. Alle Ereignisse auf diesen Strahlenwerden vom Beobachter gleichzeitig empfangen. Die Oberflache S ist im eindimensinalen Fallnur ein Punkt, dessen Bahn durch die Funktion f = 0 gegeben ist. Im eindimensionalen Fallbei Unterschallgeschwindigkeit existiert fur die Beobachtungszeit t nur einen Quellzeitpunktτ = τ ?, gegeben durch den Schnittpunkt von f und g.

Fur die Integration im mitbewegten Bezugssystem wird der Fall betrachtet, daß die durchf = 0 beschriebene Oberflache starr ist. Es liegt dann lediglich eine Translation und/oderRotation der Flache S im Raum vor. In diesem Fall kann man die Integration vereinfachen,indem man ein korperfestes Koordinatensystem ~η einfuhrt. In diesem Koordinatensystem istf nicht mehr von der Zeit abhangig. Das neue Koordinatensystem wird so definiert, daß die~η-Koordinaten mit den ~y-Koordinaten zum Zeitpunkt τ = 0 ubereinstimmen. Die beidenKoordinaten sind durch

~y = ~η + ~xs(~η, τ) (3.69)


Beobachter

Schnittpunkt

x y

τ

τ ∗

t

f

g

Abbildung 6: Bestimmung der retardierten Zeiten mit f und g

mit

~xs(~η, τ) =

∫ τ

0

~V (~η, τ ′)dτ ′

festgelegt.Nun werden die Integrale in der Losung (3.67) vertauscht und ~y mit ~η substituiert.

d3~y = Jd3~η mit J = det

(∂yi

∂ηj

)(3.70)

Dabei ist J ist die Determinante der Funktionalmatrix. Da die beiden Bezugssysteme durchreine Translation und Rotation ineinander ubergehen, ist J = 1. Ist die Substitution abge-schlossen, wird die Reihenfolge der Integrale wieder zuruckgetauscht. Fur die Modellquellefolgt:

p′(~x, t) =1

4π

∫

R3

+∞∫

−∞

1

rq(~η, τ)δ(g)dτd3~η (3.71)

mit

r(~x, ~η, τ) = | ~x− ~η − ~xs(~η, τ) | (3.72)

g(~x, t, ~η, τ) = t− r

c− τ. (3.73)

Zur Losung des inneren Integrals wird der Beobachtungsort ~x, die Beobachtungszeit t unddie Koordinate ~η festgehalten. Die Funktion g ist dann nur noch abhangig von τ . Fur einebeliebige Funktion h(τ) gilt die Regel


+∞∫

−∞

h(τ)δ(g(τ))dτ =N∑

n=1

h(τ ?n)

| dgdτ| (3.74)

mit

dg

dτ(~η, τ) = −1

c

dr

dτ|~η −1 = Mr(~η, τ)− 1. (3.75)

Dabei ist Mr die Beobachtungsmachzahl fur den Punkt mit der Koordinate ~η . Sie gibt an,mit welcher Machzahl sich der Punkt auf den Beobachter zu bewegt. τ ?

n sind die Nullstellender Funktion g. Die Summe besteht bei Unterschallgeschwindigkeit aus nur einem Wert. Mith(τ) = q/r ergibt sich fur die Losung

p′(~x, t) =1

4π

∫

R3

[q

r | 1−Mr |]

τ=τ?

d3~η. (3.76)

Die eckigen Klammern [·]τ=τ? bedeuten, daß der darin enthaltene Ausdruck zur retardiertenZeit ausgewertet werden muß. Existieren mehrere retardierte Zeiten, so mussen die Termeaddiert werden. Zur Bestimmung der retardierten Zeit τ ? wird die folgende Gleichung genutzt

c · (t− τ ?(~η)) =| ~x− ~η − ~xs(~η, τ ?(~η)) | . (3.77)

Wird jetzt wieder die flachenhafte Quellverteilung betrachtet, so muß der Betrag des Gra-dienten von f auch im ~η-Koordinatensystem gebildet werden. Da die beiden Bezugssystemedie gleiche Skalierung besitzen, sind die Betrage identisch.

| grad~yf |=| grad~ηf | (3.78)

Die Funktion f ist in dem ~η-Koordinatensystem nicht von der Zeit τ abhangig. Daher konnendie Ausdrucke mit f aus der eckigen Klammer herausgezogen werde. Zusammen mit Glei-chung (3.78) ergibt sich:

p′(~x, t) =1

4π

∫

R3

[Q

r | 1−Mr |]

τ=τ?

| grad~ηf | δ(f)d3~η (3.79)

Das Volumenintegral wird mit folgender Beziehung in ein Oberflachenintegral umgeformt:

∫

R3

h(~η)δ(f(~η))d3~η =

∫

S

h(~η)

| grad~ηf |dS(~η) (3.80)

Mit

h(~η) =

[Q

r | 1−Mr |]

τ=τ?

| grad~ηf | (3.81)

ergibt sich fur die Losung


p′(~x, t) =1

4π

∫

S

[Q

r | 1−Mr |]

τ=τ?

dS(~η). (3.82)

Dies ist die Losung als Integral uber die Oberflache in dem mitbewegten Bezugssystem.Angewandt auf die Terme der erweiterten Lighthill-Gleichung, ergibt sich die Gleichung vonFfowcs-Williams und Hawkings:

4πc2 {%′H(f)} (~x, t) =∂2

∂xixj

∫

R3

[TijH(f)

r | 1−Mr |]

τ=τ?

d3~η

+∂

∂t

∫

S

[%(ui − vi) + %0vi

r | 1−Mr | ni

]

τ=τ?

dS(~η) (3.83)

− ∂

∂xi

∫

S

[%ui(uj − vj) + Pij

r | 1−Mr | nj

]

τ=τ?

dS(~η)

Die Eigenschaften und Bedingungen fur die Gleichung von Ffowcs-Williams und Hawkingslauten:

• Prinzipiell gilt die FW-H-Gleichung in jedem nicht beschleunigt bewegten Bezugssy-stem, d.h. es muß ein Inertialsystem sein. In diesem Inertialsystem werden alle in derGleichung vorkommenden Großen gemessen.

• Alle Integrationen werden im korperfesten ~η-System durchgefuhrt. Es kann dabei aucheine beschleunigte Bewegung, z.B. Rotation, ausfuhren.

• Im Falle eines beschleunigt bewegten ~η-Systems kann dieses nicht als Bezugssystem furdie FW-H-Gleichung genutzt werden.

• Die Klammern [·]τ=τ? bedeuten, daß der Ausdruck zur retardierten Zeit τ ? ausgewertetwird. Existieren mehrere τ ?, so werden die Werte addiert.

• Zur Bestimmung der retadierten Zeiten wird Gleichung (3.77) genutzt:

c · (t− τ ?(~η)) =| ~x− ~η − ~xs(~η, τ ?(~η)) |

• Die Gleichung ist ohne Naherungen aus den nichtlinearen Erhaltungsgleichungen her-geleitet worden!


Gleichung (3.83) wird implementiert zur Berechnung des Schalles im Fernfeld. Die Losungan einem Beobachterort ~x, zu einem Zeitpunkt t, wird durch die Auswertung der zwei Ober-flachenintegrale und des Raumintegrals berechnet. Die Quellen werden zu den jeweiligenretardierten Zeiten berechnet und durch die Integration summiert. Die Auswertung einesbestimmten Zeitintervalls der instationaren Stromungsdaten ermoglicht die Berechnung desZeitverlaufs des Schalldruckes beim Beobachter und somit eine Frequenzanalyse.

3.4.1 Der FWH-Algorithmus auf einer Korperoberflache

Fur den Spezialfall, daß die Integrationsflache, beschrieben durch die Funktion f , mit derfesten Oberflache eines Korpers ubereinstimmt, sind weitere Vereinfachungen der Quelltermemoglich. Betrachtet werden wiederum die Quellterme Qi und Lij, die auf der OberflacheWerte liefern:

q(~x, t) =∂

∂t

{(%(ui − vi)ni + %0vini) | gradf | δ(f)

}

− ∂

∂xi

{(%ui(uj − vj)nj + Pijnj) | gradf | δ(f)

}(3.84)

Bei einer festen Oberflache, die sich mit der Normalengeschwindigkeit ~v bewegt, ist dieKomponente der Fluidgeschwindigkeit ~u in Richtung des Korpers gleich null. Damit ist dieGeschwindigkeit im mitbewegten Bezugssystem (ui − vi) parallel zur Korperoberflache unddamit senkrecht zum Nomalenvektor ~n. Da das Skalarprodukt zweier senkrecht aufeinander-stehenden Vektoren null ist, gilt

(ui − vi)ni = 0.

Werden noch die folgenden Abkurzungen

vn = vini

Li = Pijnj = (p− p0)ni − τijnj

eingefuhrt, so ergeben sich die vereinfachten Quellterme

q(~x, t) = +∂

∂t

{%0vn | gradf | δ(f)

}− ∂

∂xi

{Li | gradf | δ(f)

}. (3.85)

Die Oberflachenintegrale der Losung vereinfachen sich ebenso. Wird mit der linearisiertenDruck-Dichte-Beziehung p′ = c2ρ′ die Dichte mit den Schalldruck ersetzt, und die Teillosun-gen mit p′T und p′L bezeichnet, so ergibt sich:


4πp′T (~x, t) = +∂

∂t

∫

S

[%0vn

r | 1−Mr |]

τ=τ?

dS(~η) (3.86)

4πp′L(~x, t) = − ∂

∂xi

∫

S

[Li

r | 1−Mr |]

τ=τ?

dS(~η) (3.87)

Jetzt konnen die Quellterme in (3.85) auch physikalisch sinnvoll gedeutet werden. Sie re-prasentieren bewegte, flachenhafte Massen- und Impulsquellen. Die Quellterme besitzen diegleiche Form wie die Modellquellen fur bewegte, punktformige Massen- bzw. Impulsquellen(3.30 und 3.31), sind jedoch durch den Faktor | gradf | δ(f) auf der gesamten Oberflachedefiniert.

Die Große ρ0vn entspricht einer verdrangten Masse pro Flache und Zeit. Bewegt sich bei-spielsweise ein Flugel durch die Luft, so verdrangt er erst die Luft zur Seite, und danachwird sie hinter dem Flugel wieder zusammenfließen. Da jedes Oberflachenelement zu un-terschiedlichen retardierten Zeiten ausgewertet wird, liefert das Integral in p′T einen Wertungleich null. Durch die zeitliche Schwankung der verdrangten Masse entsteht Schall. Ist derFlugel symmetrisch und betrachtet man zwei gegenuberliegende Flachenelemente A1 undA2, so sind die Werte (ρ0vn)(Ai) fur die Ai vom Betrag her gleich. Durch die entgegenge-setzte Orientierung der Flachen mußten sich die Terme bei der Integration aufheben. DieWerte (ρ0vn)(Ai) werden jedoch zu unterschiedlichen retardierten Zeiten ausgewertet, da derAbstand der beiden Flachenelemente zum Beobachter unterschiedlich groß ist. Man kanndies als eine zusatzliche Phasenverschiebung, die gerade dem Langenunterschied der beidenQuellen in Richtung Beobachter entspricht, interpretieren.

Besitzt der umstromte Korper eine gewisse Dicke, so wird die Teillosung p′T einen Wertfur den Schalldruck liefern. Betrachtet man eine dunne Platte, so ist der Unterschied derretardierten Zeiten fur zwei gegenuberliegende Oberflachenelemente vernachlassigbar klein.Eine dunne Platte wird fur das Integral aus Gleichung (3.86) keinen Wert liefern. Aus diesemGrund wird dieser Anteil der Losung als “thickness-noise” bezeichnet.

Die zweite Teillosung p′L wird auch als “loading-noise” bezeichnet. Der Vektor Li entsprichteiner lokalen Kraft pro Flache, die von dem Korper auf das Fluid ausgeubt wird. Durch dieseKraft wird ein Impuls auf Medium ubertragen. Nach Newton ist “actio gleich reactio” unddie Summe der Impulse entspricht gleichzeitig dem Auftrieb bzw. Widerstand eines Flugels(in Pij sind Druck- und Scherkrafte vorhanden). p′L laßt sich somit nicht ohne weiteresverringern, es sei denn durch gleichzeitige Verringerung des Auftriebs. Es ist jedoch moglichdie Schwankung von Li zu verringern und damit den “loading-noise” p′L zu optimieren.

Als dritte Teillosung p′LH ist noch das Raumintegral mit dem Lighthillschen Spannungstensorzu betrachten.

26 3.5 Zweidimensionale Schallquellen

4πp′LH(~x, t) =∂2

∂xixj

∫

R3

[%uiuj − τij + δij(p

′ − c2%′)r | 1−Mr |

]

τ=τ?

d3~η (3.88)

Die zeitlich gemittelten Werte der Großen %uiuj werden in der Stromungslehre als Reynolds-spannungen bezeichnet. Ihr Betrag stellt eine charakteristische Große zur Beschreibung derTurbulenz einer Stromung dar. Nur in stark turbulenten Stromungen wird %uiuj einen nichtmehr zu vernachlassigen Wert zu der Gesamtlosung des Schalldrucks beisteuern. Bei derUmstromung von Korpern sind nur in einem begrenzten Gebiet VQ in unmittelbarer Nahedes Korpers die turbulenten Schwankungen groß. Im Außenbereich der Stromung wird durchden Korper im allgemeinen keine Turbulenz erzeugt. Damit kann man das Integral uber R3

auf das begrenzte Volumen VQ beschranken.

Fur den Spannungstensor τij wurden bisher noch keine Annahmen bei der Herleitung getrof-fen. Da die Stromungssimulationen auf den Navier-Stokes-Gleichungen basieren, bietet sichfur die Definition von τij in der Akustik ebenso der Schubspannungsansatz von Newton an.

τij = η

(∂uj

∂xi

+∂ui

∂xj

)(3.89)

Hier bezeichnet η die Viskositat des Fluids. Die Schubspannungen konnen bei turbulentenScherschichten große Werte annehmen. Dies ist z.B. im Freistrahl hinter einem Flugzeugtrieb-werk der Fall. Eine Verringerung des Larms wird durch Reduzierung der Austrittsgeschwin-digkeit des Strahls erreicht. Bei der Schallentstehung eines Flugzeugtriebwerk dominiert dieTeillosung p′LH gegenuber p′T und p′L. Im Allgemeinen ist der Lighthillterm kleiner als der“loading”- oder “thickness”-Noise.

Der Term δij(p′ − c2%′) ist nur bei großen Temperaturschwankungen von Bedeutung. In

diesem Fall gilt die linearisierte Druck-Dichte-Beziehung nicht mehr.

3.5 Zweidimensionale Schallquellen

Bisher wurden nur dreidimensionale Schallquellen behandelt. Losungen fur zweidimensionaleSchallquellen, d.h. konstante Quellstarke in x3-Richtung, sind aus den dreidimensionalenLosungen herleitbar.

3.5.1 Die zweidimensionale Greensche Funktion

Mit einer zweidimensionalen Greenschen Funktion ist die Losung fur zweidimensionale Schall-quellen als Integral uber die Schallquelle multipliziert mit der Greenschen Funktion darstell-bar.

Die inhomogene Wellengleichung mit dem Quellterm q auf der rechten Seite der DGL lautet

(1

c2

∂2

∂t2−∆

)φ = 4πq(~x, t). (3.90)

Die Losung dargestellt mit Hilfe der Greenschen Funktion aus Kapitel 3.2.1 lautet

3.5 Zweidimensionale Schallquellen 27

φ(~x, t) =

+∞∫

−∞

∫

R3

G(~x, t, ~y, τ)q(~y, τ)d3~ydτ (3.91)

mit

G(~x, t, ~y, τ) =δ(t− r

c− τ)

4πrmit r =| ~x− ~y | . (3.92)

Die zweidimensionale Quellverteilung Q ist abhangig von x1, x2, t und konstant in x3-Richtung.

q = Q(x1, x2, t) (3.93)

Die zweidimensionale Losung lautet dann

φ(x1, x2, t) =

+∞∫

−∞

+∞∫

−∞

+∞∫

−∞

G2D(x1, x2, t, y1, y2, τ)Q(y1, y2, τ)dy1dy2dτ. (3.94)

Die neue Greensche Funktion G2D ist nicht mehr von x3 und y3 abhangig und ist durch dasIntegral

G2D =

+∞∫

−∞

G(~x, t, ~y, τ)dy3 =

+∞∫

−∞

δ(t− rc− τ)

4πrdy3 (3.95)

gegeben. Durch Substitution des Argumentes der δ-Funktion mit σ = t − rc− τ folgt die

Beziehung zwischen den Differentialen

dy3 =cr

x3 − y3

dσ. (3.96)

σ wachst dabei von −∞ bis zum Maximum bei t− Rc−τ und fallt danach wieder bis −∞ ab.

Das Maximum wird erreicht wenn sich die Quelle und der Beobachter auf gleicher Hohe inx3-Richtung befinden. Wie in Abbildung 7 gezeigt, ist in der x1-x2-Ebene des Beobachters derAbstand von Beobachter zur Quelle in dieser Ebene am geringsten. Der im zweidimensionalenRaum gebildete Abstand wird mit R bezeichnet. Mit zunehmenden x3, in positiver wie auchnegativer Richtung, wird der Abstand r =

√R2 + s2 großer, analog dazu σ kleiner.

Die Symmetrie des Problems ermoglicht die Teilung des Integrals, wobei nach Vertauschender Integrationsgrenzen des einen Integrals und dem Einfuhren des Betrages von x3−y3 zweiidentische Teilintegrale entstehen. Es ergibt sich

G2d =c

2π

t−Rc−τ∫

−∞

δ(σ)

| x3 − y3 |dσ. (3.97)

Falls die obere Integrationsgrenze t− Rc− τ kleiner null ist, dann ist der Wert des Integrals

gleich Null. Andernfalls ergibt sich der Wert des Quotienten an der Stelle σ = 0. Damiterhalt man die Greensche Funktion fur den zweidimensionalen Fall


R

r

r

s

x 1

x 2

x 3

Beobachter

zweidimensionale diskrete Punktquelle im R

x 1

x 2

R

3

2D:

3D:

zweidimensionale Punktquelle im R 2

Abbildung 7: Zweidimensionale diskrete Quellverteilung

G2D =

{c

2π√

c2(t−τ)2−R2fur t− R

c− τ > 0

0 sonst.(3.98)

Die Losung fur zweidimensionale Schallquellen ist damit wieder als Integral uber die Quell-terme und der neuen Greenschen Funktion darstellbar.Definiert man als Quellterm eine einfache zweidimensionale Punktquelle, die einen δ-Puls imKoordinatenursprung aussendet, so erhalt man als Losung die zweidimensionale GreenscheFunktion G2D. Die Losung ist in Abbildung 8 dargestellt und lautet

φ(x1, x2, t, 0, 0, 0) =c

2π√

c2t2 −R2fur t− R

c> 0. (3.99)

G2D ist analog zum dreidimensionalen Fall als Elementarwelle fur 2D-Quellen interpretierbar.Die Losungen des zwei- bzw dreidimensionalen Fall unterscheiden sich jedoch stark. Wahrendim dreidimensionalen Fall der δ-Puls an jedem Beobachtungsort zur retardierten Zeit wiedereinen δ-Puls bewirkt, bleibt die Signalform im zweidimensionalen Fall nicht erhalten. Zurretardierten Zeit hat das Signal eine Singularitat, da es dort nicht definiert ist. Fur jedenZeitpunkt nach der retardierten Zeit registriert der Beobachter im Abstand R ein Signal vonder Quelle. Der Signalverlauf ist monoton fallend und geht erst fur t →∞ gegen null.

Ein unendlich kurzer 2D-Puls erzeugt ein unendlich lang abklingendes Signal.

Das Signal nimmt mit zunehmendem Abstand mit 1√R

ab. Linienhafte zweidimensionaleSchallquellen, z.B. eine stark befahrene Autobahn, sind in großer Entfernung lauter als einegleichstarke Punktquelle, da im 3D-Fall die Signalstarke mit 1

rabnimmt.

Zur effektiven Berechnung einer beliebigen zweidimensionalen Quellverteilung ist die einfachezweidimensionale Greensche Funktion aus Gleichung (3.98) jedoch nicht geeignet. Es ist nichtpraktikabel, alle Quellen fur alle Zeiten vor einem bestimmeten Zeitpunkt zu bestimmen.


0 ct1

ct2

ct3 R

φSignal zur Zeit t = t

1

Signal zur Zeit t = t2

Signal zur Zeit t = t3

1/R-0.5

2D-δ-Puls im Ursprung zur Zeit t = 0

Abbildung 8: Zweidimensionale Elementarwelle zu drei unterschiedlichen Zeiten ti

3.5.2 Die harmonische zweidimensionale Schallquelle

Betrachtet man eine harmonische zweidimensionale Schallquelle Q im Koordinatenursprungmit der Zeitfunktion f

Q(x1, x2, t) = δ(x1)δ(x2)f(t) mit f(t) = Aeiωt (3.100)

so ergibt sich die Losung

φ(R, t) =Ac

2π

t−Rc∫

−∞

eiωτ

√c2(t− τ)2 −R2

dτ. (3.101)

Das Integral vereinfacht sich mit folgender Substitution

ξ =c(t− τ)

R(3.102)

mit den Beziehungen fur die Aufspaltung des Exponentialterms und des Differentials

τ = t− Rξ

cund dτ = −R

cdξ. (3.103)

Es ergibt sich

φ(R, t) =Aeiωt

2π

∞∫

1

e−iωξ Rc√

ξ2 − 1dξ. (3.104)

Die Zeitabhangigkeit steht vor dem Integral, d.h. die zweidimensionale harmonische Quellefuhrt zu einem harmonischen Zeitsignal, skaliert mit dem Wert des Integrals. Das Integralbesitzt die gleiche Form wie die Hankel-Funktion nullter Ordnung zweiter Art


H(2)0 (s) =

2i

π

∞∫

1

e−isξ

√ξ2 − 1

dξ (3.105)

= J0(s)− iY0(s).

Der Funktionswert von H(2)0 (s) ist komplex. Der Realteil entspricht der Besselfunktion nullter

Ordnung J0, und der Imaginarteil entspricht der Neumannfunktion Y0. Die Hankel-Funktionist nicht analytisch bestimmbar, sie kann jedoch mit Hilfe der Bessel- und Neumannfunktiondurch Reihen dargestellt werden. Damit sind die Werte numerisch mit der gewunschterGenauigkeit berechenbar. Die Funktionen sind sind in Abbildung 9 dargestellt. Sie besitzeneinen wellenformigen Verlauf. Die Hankelfunktion beschreibt eine nach außen laufende Welleeiner zweidimensionalen Quelle. Die einhullende Kurve ist eine Funktion proportional zu 1√

s.

0 10 20 30 40s

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1J

0(s), Y

0(s)

J0(s)

Y0(s)

(-)0.8/s0.5

Abbildung 9: Bessel- und Neumannfunktion nullter Ordnung

Die Losung fur nach außen laufende Wellen fur zweidimensionale Schallquellen lautet dann

φ(R, t) =Aπ

2iceiωtH

(2)0

(ωR

c

). (3.106)

Existieren mehrere Quellen, so ist die Gesamtlosung durch Uberlagerung der Einzellosungennach Gleichung (3.106) bestimmbar.

3.5.3 Zweidimensionale Greensche Funktion in der Frequenzebene

David P. Lockard [8] stellt noch eine weitere Variante der zweidimensionalen GreenschenFunktion vor. Er formuliert die Wellengleichung in einem bewegten Koordinatensystem undfuhrt dann eine Fourier-Transformation der Wellengleichung durch. In Lockards Algorithmus


werden dabei fur jeden einzelnen Quellpunkt die Quellen fur einen festen Zeitraum ausge-wertet und danach die Fouriertransformationen mit Hilfe einer FFT (Fast Fourier Trans-formation) berechnet. Zwingende Voraussetzung ist dabei der konstante Abstand zwischenBeobachter und Quellen. Die dazu passende Greensche Funktion lautet

G(x, y, ξ, η) =i

4βe

Mkxβ2 H

(2)0

( k

β2

√x2 + β2y2

)(3.107)

mit ξ und η als Quellkoordinaten, der Beobachterposition in x und y und den folgendenBeziehungen

x = (x− ξ) cos θ + (y − η) sin θ

y = −(x− ξ) sin θ + (y − η) cos θ

tan θ =V

U

β =√

1−M2

M =

√U2 + V 2

c

k =ω

cω = 2πf.

Die Pegel fur die einzelnen Frequenzen (bzw den Wellenzahlen k) fur jeden Quellpunktwerden mit Hilfe der speziellen Greenschen Funktion aufsummiert.Die Bestimmung der Losung ist also in der Zeit- oder Frequenzebene moglich.

32 4.1 Zweidimensionale Implementierung des FW-H

4 Numerisches Verfahren

4.1 Zweidimensionale Implementierung des FW-H

Die Gleichung von FW-H kann in der Frequenz- und Zeitebene angegeben werden.

Fur die Implementierung in der Frequenzebene mussen alle Terme der FW-H-Gleichungumgeformt werden, wobei zusatzliche Einschrankungen gemacht werden. So wird bei demvon Lockard [8] vorgestellten Algorithmus ein konstanter Abstand von allen Quellpunk-ten zum Beobachter vorausgesetzt, und das Bezugssystem fur den FW-H ist mitbewegt.Das FW-H-Bezugssystem ist bei Lockard identisch mit dem Koordinatensystem, in dem dieKontrollflache f definiert ist und damit auch die Integrale uber die Quellen gebildet wer-den. Die Quelle kann somit nur rein translatorische, gleichformige Bewegungen ausfuhren.Dies schrankt die moglichen Anwendungsfalle ein. Da in den Herleitungen der GreenschenFunktion nach Lockard die zweidimensionale Quellverteilung eingeflossen ist, sind die vor-hergesagten Pegel im Vergleich zu realen Messungen immer zu hoch. Der Algorithmus vonLockard wird aus diesen Grunden hier nicht verwendet.

Da in der Realitat die umstromten Korper, z.B. der Flugel eines Flugzeuges, in der drittenRaumrichtung endlich sind, der zweidimensionale FW-H die Quellen jedoch in die Unend-lichkeit fortsetzt, wird der Pegel immer zu hoch vorhergesagt. In dieser Arbeit soll jedochein Algorithmus fur die Schallvorhersage von realen, umstromten Korpern auf Basis von2D-CFD-Daten implementiert werden. Die Verwendung der zweidimensionalen GreenschenFunktion, wie in Kapitel 3.5 beschrieben, ist dann jedoch nicht moglich.

Eine mogliche Losung ist, nur begrenzt in x3-Richtung zu rechnen. Die Berechnung desSchalldruckes erfolgt quasi-zweidimensional, indem die 2D-Daten aus der Stromungsimula-tion bis zu einem festen Wert in x3-Richtung identisch fortgesetzt werden. In diesem Fallkann der dreidimensionale FW-H mit der entsprechenden, einfachen Greenschen Funktionverwendet werden. Dazu muß fur die Akustik eine geignete Diskretisierung in x3-Richtunggewahlt werden. Die Situation ist in Abbildung 10 dargestellt.

Aus der zweidimensionalen Stromungssimulation erhalt man entlang einer Linie an diskretenPunkten i Werte fur den Druck und die Geschwindigkeit. Mit der Elementbreite ∆z in x3-Richtung entstehen die Flachenelemente Ai. Der Mittelpunkt dieser Flachenelemente liegtbei x3 = 0. Die Position des Beobachters in x3-Richtung ist ebenfalls bei x3 = 0. Da ∆z in dergleichen Großenordnung wie die Diskretisierung der Stomungssimulation liegen sollte, werdenKorper mit einer bestimmten Ausdehung durch Aneinanderreihen identischer Elemente in x3-Richtung approximiert. Die Mittelpunkte der Elemente sind durch die Formel x3(j) = j ·∆zmit j = 0, . . . , nz gegeben. Wird beispielsweise ∆z = 0, 02m gewahlt und die Anzahl derElemente nz in x3-Richtung mit 25 festgelegt, so ist die Ausdehnung des Korpers max(x3) =nz ·∆z + ∆z

2= 0.51m. Wird im Algorithmus noch die Symmetrie das Problems ausgenutzt,

so verdoppelt sich die Gesamtausdehnung lz des Korpers in x3-Richtung auf 1.02m.

lz = 2∆z(nz +1

2) (4.1)

Der FW-H nach Gleichung (3.83) ist in der Zeitebene formuliert. Es ist moglich, den FW-Hin zwei verschiedenen Varianten zu implementieren.

4.1 Zweidimensionale Implementierung des FW-H 33

x1

x3

x2

z∆

l = 2 z (nz+1/2) ∆z

i

x1i

x2i

hier nz = 3

k

k

x3 x3k = nz z ∆.

x2k x2

i

x1k x1

i

=

=

.

p = pi k

i, k − Element i bzw k

Abbildung 10: Fortfuhren der CFD-Daten in x3 −Richtung fur die 2D-Akustik

1. Die Quellterme fur alle Punkte werden fur eine feste Quellzeit berechnet und die re-tardierten Zeiten aller Quellpunkte bestimmt. Danach ist es moglich, die Quelltermezum Schalldruck zu den richtigen Beobachtungszeiten aufzusummieren.

2. Fur eine feste Beobachtungszeit werden die retardierten Zeiten aller Elemente be-stimmt. Die Quellterme werden dann fur die korrekten Quellzeiten berechnet.

Nur Variante zwei ist konsistent zu implementieren.Die Numerik rechnet immer diskret. Die Beobachterzeit wird mit der gleichen Auflosungwie die Quellzeit diskretisiert ∆t = ∆τ . Die in Variante 1 bestimmten Beobachtungszeiten,stimmen im allgemeinen nicht mit den diskreten Zeitpunkten ti uberein (∆t = ti − ti−1).Die errechneten Werte fur den Schalldruck der einzelnen Elemente mussen interpoliert wer-den. Bei der ersten Variante ist jedoch bei unabhangig voneinander bewegtem Beobachterund Korper eine lineare Interpolation auf die dieskreten Beobachterzeiten nicht moglich.Zusatzlich besteht die Moglichkeit, daß die berechneten Beobachtungszeiten gehauft um denZeitpunkt ti−ti−1

2auftreten. Werden die Quellterme ungleichmaßig auf ti und ti−1 aufgeteilt

und interpoliert, so werden die Schalldrucke zu diesen Zeitpunkten zu niedrig bzw. zu hochvorhergesagt.Bei Variante 2 funktioniert die lineare Interpolation der Quellwerte auf Basis der Werte vondrei Quellzeiten (Fehler max O(h2)) auch bei unabhangiger Bewegung von Beobachter undKorper, und es wird garantiert, daß fur jeden Beobachterzeitpunkt von jedem Element nurein Quellwert berucksichtigt wird.

Die zweidimensionale Implementierung des FW-H basiert auf einem vorhandenen 3D-FW-H-Code von D. Eschricht [11], entwickelt am HFI der TUB. Der Code nutzt die beobach-terzeitbasierte Variante 2 des FW-H-Algorithmus. Der Code ist ausfuhrlich getestet und

34 4.2 Schnittstelle zwischen Stromungssimulation und FW-H

validiert.Bei Nutzung des dreidimensionalen FW-H auf der Basis von Daten aus zweidimensionalenStromungssimulationen mussen im kompletten Programmcode die Schnittstellen nach Außenund innerhalb des Programms und die Routinen zur Berechnung der Oberflachendaten mo-difiziert werden. Die Berechnung der Quellterme und die Bestimmung der Quellzeiten sindidentisch zum 3D-Code. Die Werte aus der zweidimensionalen Stromungssimulation werdenin dem Programm den Elementen im R3 zugeordnet.

4.2 Schnittstelle zwischen Stromungssimulation und FW-H

Wird die Umstromung von Korpern simuliert, erhalt man im gesamten diskretisierten Re-chengebiet die instationaren Felder der Fluidgeschwindigkeit ui, des Druckes p und weitereinteressierende Stromungsgroßen, z.B. die turbulente kinetische Energie. Zur Berechnungdes Schalldruckes im Fernfeld mit dem FW-H nach Gleichung (3.83) werden fur die zweiTeillosungen pL und pT nur Werte auf einer geschlossenen Kontrollflache S, beschriebendurch die Funktion f , um den Korper benotigt. Dabei kann die Kontrollflache auch identischmit der Korperoberflache sein. Alle akustischen Quellen innerhalb und auf der Kontrollfachewerden damit erfasst. Fur den Ligthillterm pLh der Losung muß das Integral der Quelltermeim gesamten Raum R bekannt sein und auch ausgewertet werden. Fur den gesamten Raumkann die Stromungssimulation jedoch nicht die Werte liefern.Fur die Berechnung des Schalldruckes mit dem FW-H mussen im gesamten Gebiet die Schall-drucke schon bekannt sein, auch in der unmittelbaren Umgebung des Beobachters. Einesinnvolle Anwendung der Losung nach Ffowcs-Williams und Hawkings ist so nicht moglich.

Bei der Betrachtung des Stromungsfeldes um Korper ist noch eine weitere Interpretationder Losung moglich: Bewegt sich ein Korper, beipielsweise ein Flugzeug, durch die Luft, sowird es nur in seiner unmittelbaren Umgebung das ruhende Medium storen, d. h. zeitlicheAnderungen von Geschwindigkeiten, Druck usw. verursachen. Mit zunehmender Entfernungwird der Einfluß des Flugzeuges auf die Umgebung immer geringer. Die akustisch wirksamenQuellen sind somit auf ein endliches Quellvolumen VQ beschrankt. In diesem Volumen be-findet sich der umstromte Korper. Befinden sich alle Quellen in dem Volumen, so kann dasRaumintegral uber R in ein Integral des begrenzten Volumen VQ umgewandelt werden. DasRaumintegral ist dabei nur außerhalb der Kontrollflache S definiert. Es mussen so fur jedenZeitschritt nicht die kompletten Felder aus der Stromungssimulation gespeichert werden,sondern nur die Werte auf der Oberflache S und in dem Volumen VQ.

Je nach Lage der Kontrollflache und Große des Quellvolumens gibt es vier mogliche Vari-anten, verdeutlicht in Abbildung 11, die relevanten Stromungsdaten fur die Akustikberech-nungen zu extrahieren.

1. Die Kontrolloberflache S ≡ S1 ist identisch mit der Korperoberflache. Der pLh-Termwird vernachlassigt, d.h. es wird kein Volumenintegral ausgewertet.

Der Fehler des Verfahrens ist groß, da alle Effekte in der Stromung, z.B interagieren-de Wirbel oder starke Scherschichten, vernachlassigt werden. Hauptargument fur die

4.2 Schnittstelle zwischen Stromungssimulation und FW-H 35

CFD − Rechengebiet

VQ

u0

SS S1 23

V − akustisch wirksames GebietQ

S (f) − moegliche Kontroflaechen definiert durch fi i

Abbildung 11: Wahl von Kontrollflachen und Lage des Quellvolumens

Beschrankung auf die Korperoberflache ist, daß bei den meisten Stromungskonfigura-tionen die Quellterme auf der Oberflache großer als die Quellterme in der Stromungsind. Vorteil sind die vereinfachten Oberflachenintegrale (siehe Gleichungen (3.86) und(3.87)), der numerische Aufwand sinkt. Ebenso mussen nur noch die Werte fur denDruck aus der Stromungssimulation auf der Korperoberflache gespeichert und ausge-wertet werden.

Diese Variante kommt sehr haufig zur Anwendung, da bei komplizierten, rotierendenProblemen die Bestimmung der retardierten Zeiten aufwandig ist. Jeder eingesparteQuellpunkt verringert den numerischen Aufwand der Gesamtrechnung.

2. Die Kontrolloberflache S ≡ S1 ist identisch mit der Korperoberflache. Der pLh-Termwird berechnet, d.h. es wird das Volumenintegral ausgewertet.

Wird das Volumen ausreichend groß gewahlt, werden alle akustisch wirksamen Quellenerfasst. Der Fehler des Algorithmus geht gegen null. Diese Variante wird sehr seltenangewandt, da die Anzahl der Quellpunkte stark anwachst und damit die Rechenzei-ten lang werden (mindestens eine Großenordnung mehr gegenuber Variante 1). DieImplementierung des kompletten Lighthillterms ist mit erheblichen Mehraufwand ver-bunden, da zweifache Ortsableitungen der Quellterme gebildet und integriert werdenmussen. Die Erfassung der instationaren Werte in einem Raumgebiet vergroßert dieDatenmengen im Vergleich zu Variante 1, 3 und 4 stark.

Die Anwendung der Variante mit dem Lighthillterm ist selten und wird dann meistnur zu Vergleichszwecken herangezogen. Die Genauigkeit ist gegenuber den anderenVarianten am hochsten ist.

3. Die Kontrolloberflache S ≡ S2 schließt das Quellvolumen VQ ein. Es wird kein Raum-integral ausgewertet.

Da das Gebiet mit den aktiven akustischen Quellen von S2 eingeschlossen ist, werdendurch Auswertung der Quellen auf der Kontrolloberflache auch alle Quellen innerhalb


von S2 erfasst. Der Fehler des Algorithmus geht analog zu Variante 2 gegen null. DerStromungsloser muß jedoch in der Lage sein die Schallausbreitung der realen Quellenauf der Korperoberflache und im Volumen VQ bis zur Kontrolloberflache korrekt zu be-rechnen. Stromungsloser weisen im allgemeinen eine hohe numerische Dissipation auf.Wird die Kontrolloberflache weit entfernt von dem umstromten Korper gelegt, wird einGroßteil der kleinen Druckschwankungen - die in Wandnahe verhanden sind - dissipiertsein und damit gar nicht mehr erfaßt. Die Ausgangsdaten fur die Akustikrechnung sinddamit schon fehlerbehaftet. Der Fehler pflanzt sich dann in dem FW-H Algorithmusfort, die Ergebnisse stimmen nicht mehr.

4. Die Kontrolloberflache S ≡ S3 wird in die Nahe des Korpers gelegt. Es wird keinRaumintegral ausgewertet.

Wird die Kontrolloberflache dichter an den Korper herangelegt, so sind die Schalldruck-amplituden innerhalb von S3 großer und der Stromungsloser muß eine nicht so großeEntfernung uberbrucken. Die Diskretisierung des Rechengebietes ist in unmittelbarerNahe des Korpers feiner, die Stromungssimulation ist dann nicht so stark fehlerbehaf-tet. Der Fehler bei der Vorhersage der Quellterme auf der Kontrolloberflache verringertsich. Der Fehler durch Vernachlassigung von Quelltermen außerhalb von S3 wird je-doch großer. Die richtige Wahl der Oberflache ist ein Kompromiß und muß fur jedeStromungskonfiguration individuell festgelegt werden.

Ausfuhrliche Untersuchungen zur Wahl der Kontrolloberflache sind beispielsweise vonBrentner und Farassat [1] durchgefuhrt worden.

Bei dieser Variante werden nur die Werte auf der Oberflache S benotigt, die Anzahl derdiskreten Flachenelemente wird im allgemeinen großer als bei Variante 1 sein, da dieGesamtflache großer wird. Da die Auswertung der Quellterme bei einer durchlassigenOberflache aufwandiger als bei einer festen Korperoberflache ist, wird der numerischeAufwand im Verhaltnis zu Variante 1 steigen, jedoch in der gleichen Großenordnungliegen. Die Genauigkeit der Ergebnisse ist ahnlich gut wie bei Variante 2 und 3.

Bei einer geschickten Wahl der Kontrolloberflache vereint Variante 3 die Vorzuge deranderen Varianten bei nur moderater Steigerung des numerischen Aufwandes.

Im Rahmen dieser Arbeit ist die Kontrolloberflache entsprechend der Variante 4 frei wahlbar.Als Spezialfall ist es auch moglich nur die Quellterme auf der festen Korperoberflache zuberechnen.

4.2.1 Wahl der Koordinatensysteme beim FW-H

Die richtige Wahl der Koordinatensysteme ist entscheidend fur die effiziente Berechnung desSchalldruckes mit dem FW-H. Ausgangspunkt sind die Werte der numerischen Stromungssi-mulation im ~ζ-Koordinatensystem. In diesem Koordinatensystem ruht der Korper und wirdvon rechts mit einer vorgegebenen Geschwindigkeit u0 angestromt, wie in Abbildung 12averdeutlicht.In dem ~ζ-Koordinatensystem werden auf der Kontrolloberflache S die Geschwindigkeits-komponenten u, v und der Druck p fur alle Zeitschritte herausgeschrieben. Die Wahl desBezugskoordinatensystem fur den FW-H ist an bestimmte Bedingungen geknupft:

4.2 Schnittstelle zwischen Stromungssimulation und FW-H 37

Rechengebiet CFD

Koerper const. in

u

u, v, p (t) auf f

0

a)

b)

ζ

ζ 1

2

ζ

Schnittstelle

FW−H

v = uv 0

t = tt = t

01

x(t)η(t )0η(t )1

const. in ζ

const. in ηη

v

= t + t 0 ∆

∆ tv .

Bahn des Koerpers

f

fMedium in Ruhe

Abbildung 12: Koordinatensysteme bei a) CFD und b) FW-H

Das Bezugskoordinatensystem muß ein Inertialsystem sein, d.h. es darf sich nicht beschleu-nigt bewegen. So ware es moglich, das ~ζ-System auch fur den FW-H zu nutzen und auchdie Integrationen in diesem auszufuhren. Der Wellenoperator, wie in Kapitel 3.3 hergeleitet,gilt jedoch in diesem Bezugssystem nicht. Der Wellenoperator mußte mittels einer Galilei-Transformation fur die Schallausbreitung in bewegten Medien umformuliert werden. DieTransformation hat dann auch Auswirkung auf die Form der Quellterme. Ohne weitere Maß-nahmen ist diese Wahl des Bezugssystems nicht moglich.

Eine sinnvollere Wahl ist die in Abbildung 12b gezeigte. In dem Bezugssystem ~x fur denFW-H ruht das Medium. Der Korper bewegt sich mit v = u0 in negative x1-Richtung.Alle Geschwindigkeiten werden in dem ~x-Koordinatensystem angegeben, mussen also umu0 verringert werden. Das Koordinatensystem ~η, in dem die Integrationen gebildet werden,bewegt sich mit dem Korper in die negative x1-Richtung. Zum Zeitpunkt t = t0 sind ~x- und~η-Koordinatensystem identisch. Im ~η-System bewegen sich der Korper und die Kontrollflachef wie im ~ζ-Koordinatensystem der Stromungssimulation nicht. Die Position des Beobachtersmuß im ~x-System zur Zeit t = t0 definiert werden. Soll der Schalldruck immer in einem


festen Abstand und Richtung zum umstromten Korper ausgewertet werden, muß sich derBeobachter ebenfalls mit vB = u in negative x1-Richtung bewegen.Folgt auf diese Weise die Wahl der Koordinatenssteme, kann der FW-H-Algorithmus in derForm von Gleichung (3.83) zur Berechnung des Schalldruckes benutzt werden. Die freien Kon-stanten %0, p0 und c sind dann analog zur Herleitung mit den Ruhewerten (ohne Schall) bzw.der Schallgeschwindigkeit in Luft festgelegt. Nur unter diesen Voraussetzungen beschreibtder FW-H die Schallerzeugung und die Schallausbreitung im Fernfeld korrekt.

4.2.2 Konvertierung der Stromungsdaten fur die Akustik

Fur die Stromungssimulationen wird im Rahmen dieser Arbeit das Programm ELAN2 ver-wendet. Der Code rechnet dimensionslos, d.h. alle Strommungsgroßen sind durch entspre-chende Referenzgroßen dimensionslos gemacht. Die Navier-Stokes-Gleichung fur die Ge-schwindigkeitskomponenten ui lauten dann beispielsweise

∂ui

∂t+ uj

∂ui

∂ζj

= − ∂p

∂ζi

+1

Re

∂2ui

∂ζ2j

(4.2)

mit den Definitionen der dimensionslosen Großen

ζi =xi

ll - charakteristische Lange

ui =ui

u∞u∞ - Anstromgeschwindigkeit

t =uit

l

p =p

%u2∞

Re =u∞l

%νmit ν =

η

%.

Dabei ist ν die kinematische Viskositat.

Die Ausgabe der Daten aus dem CFD-Programm ELAN2 erfolgt zur Zeit formatiert. Fur je-den Zeitschritt wird eine Datei angelegt, in der die Koordinaten jedes Elementes und die ent-sprechenden Druck und Geschwindigkeitsdaten. Da fur die Auflosung von hohen Frequenzensehr kleine Schrittweiten, fur eine ausreichende Freqeunzauflosung ∆f eine moglichst langeauswertbare Gesamtsimulationszeit benotigt wird, werden sehr viele Augabedateien erzeugt(zwischen 10000 und 20000). Der FW-H-Algorithmus muß sehr haufig auf die Stromungs-daten zugreifen. Die Verwendung von unformatierten, binaren Eingabedateien beschleunigtden Algorithmus stark.Im Gegensatz zur Stromungssimulation ist in der Akustik die Verwendung von dimensi-onsbehafteten Werten von Vorteil, da man an einem “richtigen” Pegel interessiert ist unddimensionslose Frequenzen nicht gut interpretierbar sind.

4.3 Implementierung des FW-H in C2Noise 39

Aus diesen Grunden wurde im Rahmen dieser Arbeit das Tool mdat2bin entwickelt. Es liestdie formatierten Ausgabedateien der Stromungssimulation nacheinander ein und erzeugt dreineue Dateien fur das Programm C2Noise. Es werden eine formatierte Datei mit den Koor-dinaten der Kontrolloberflache und zwei unformatiert Dateien, jeweils fur Druck und Ge-schwindigkeiten auf der Oberflache, erzeugt. Dabei berechnet mdat2bin aus den dimensionslo-sen Stromungsdaten gleichzeitig die dimensionsbehafteten Werte. Die Datei akustik.paramenthalt die Eingabe- bzw. Ausgabeparameter der Stromungssimulation und der Akustik zurkorrekten Dimensionalisierung sowie den Testfallnamen. mdat2bin erzeugt des weiteren ei-ne Datei constants fur das Programm C2Noise mit den Dateinamen der unformatiertenGeschwindigkeits- und Druckdateien, des Koordinatenfiles, sowie den Parametern der Aku-stik.Mit dem Programm C2Noise, der Parameterdatei constants und der Datei observers, indenen die Koordinaten der Beobachter stehen, ist die Berechnung des Schalls im Fernfeldmit dem FW-H-Algorithmus moglich.

4.3 Implementierung des FW-H in C2Noise

In den folgenden Kapiteln wird gezeigt, wie der Algorithmus von FW-H im ProgrammC2Noise umgesetzt wurde. Dabei wird auch auf die Schnittstelle zwischen der Stromungssi-mulation und dem Akustikprogramm eingegangen.

4.3.1 Umwandlung der Orts- in eine Zeitableitung

Ortsableitungen sind numerisch schlecht zu implementieren. Es gibt die Moglichkeit, dieOrtsableitung in eine Zeitableitung umzuformen. In einem nicht beschleunigt bewegten Be-zugssystem ~η (analog zu Kapitel 3.4) ist die inhomogene Wellengleichung mit der flachen-haften Quellstarkeverteilung Wi gegeben:

(∂2

∂t2− c2∆

)φ =

∂

∂xi

{Wi | gradf | δ(f)

}(4.3)

Die Losung ist als Integral uber die Oberflache S gegeben:

φ(~x, t) =1

4π

∂

∂xi

∫

S

[Wi(~η, τ)

r | 1−Mr |]

τ=τ?

dS(~η) (4.4)

Dann ist die Umformung der Ortsableitung mit folgendem Zusammenhang moglich:

∂

∂xi

∫

S

[Wi

r | 1−Mr |]

τ=τ?

dS(~η)

(4.5)

= − ∂

∂t

∫

S

[Wr

cr | 1−Mr |]

τ=τ?

dS(~η)−∫

S

[Wr

r2 | 1−Mr |]

τ=τ?

dS(~η)

40 4.3 Implementierung des FW-H in C2Noise

mit

Wr = Wi(x− y)i

r. (4.6)

Bie Große Wr ist die Komponente des Vektors ~W in Richtung des Beobachters. Der Ausdruck(x−y)i

rstellt den normierten Vektor dar, der vom Quellpunkt in Richtung des Beobachters

zeigt, wie in Abbildung 13 gezeigt. Die skalare Große Wr ist je nach Winkel zwischen demVektor Wi und dem Vektor (x−y)i

runterschiedlich groß. Ist Wi beispielsweise in der gleichen

Richtung wie u orientiert, ergibt sich fur eine auf den Beobachter zubewegte Quelle einpositiver Wert. Analog gilt fur eine wegbewegte Quelle, daß der Wert Wr negativ wird.Befinden sich die Quelle und der Beobachter auf einer Hohe, so ist Wr gleich null.

Die Bewegung der Quelle verandert, je nach effektiver Geschwindigkeit zum Beobachter, denWert der Quelle, die der Beobachter “sieht”. Bei der an einem Beobachter vorbei bewegtenQuelle andert sich nicht nur die Frequenz, der sogenannte Dopplereffekt, sondern auch dieAmplitude des empfangenen Signals.

.

.

.

(x−y)r

i (t )1

v

W (t ) i 1

W (t ) i 2 W (t ) i 3

W (t ) > 0r 1

(x−y)r

i (t )2

(x−y)r

i (t )3

W (t ) = 0r 2

W (t ) < 0 r 3

Beobachter am Ort x i

Bahn der Quelle y (t)i

Abbildung 13: Die effektiven Werte bei bewegter Quelle

Die Summe uber die drei Integrale der Ortsableitungen wird mit Gleichung (4.6) in ein Inte-gral ohne Ableitung und ein zweites mit Zeitableitung umgeformt. Es ergibt sich insgesamteine Vereinfachung. Desweiteren lassen sich die Teillosungen p′T und p′L teilweise zusam-menfassen. Angewandt auf die Gleichung von FW-H (3.83) ergibt sich fur die betrachtetenOberflachenintegrale fur eine durchlassige Kontrolloberflache:


4πp′(~x, t) =∂

∂t

∫

S

[%0Un

r(1−Mr)

]

τ=τ?

dS(~η)

︸︷︷︸p′T

+∂

∂t

∫

S

[Lr

cr(1−Mr)

]

τ=τ?

dS(~η) +

∫

S

[Lr

r2(1−Mr)

]

τ=τ?

dS(~η)

︸︷︷︸p′L

(4.7)

mit den Beziehungen

un = uini

vn = vini

Un = Uini

=

(1− %

%0

)vn +

%un

%0

Lr =(x− y)i

rLi

Li = %ui(un − vn) + Pin

Pin = Pijnj mit Pij = (p− p0)δij − τij

= (p− p0)ni | τij vernachlassigt

Mr =vr

cr = | ~x− ~η − ~xs(~η, τ) |

p′ = c2%′ linearisierte Druck-Dichte-Beziehung.

Die Bildung des Betrages beim Ausdruck 1−Mr wird weggelassen, da das Programm nur furUnterschallgeschwindigkeit ausgelegt ist. In diesem Fall ist immer gewahrleistet, das 1−Mr

positiv ist.Diese Formulierung hat jedoch einen Nachteil. Um den Schalldruck an einem Ort zu einembestimmten Zeitpnkt zu berechnen, mussen neben den drei Integralen zur retardierten Zeitauch die Integrale fur die Quellzeitpunkte vor bzw. hinter der retardierten Zeit bestimmtwerden. Diese Werte werden fur die Berechnung der Zeitableitungen benotigt.Es ist moglich, die Zeitableitungen mit der folgenden Beziehung in die Integrale hineinzuzie-hen.


∂

∂t=

[1

1−Mr

∂

∂τ

]

ret

(4.8)

Die Ableitungen nach τ werden mit folgenden Beziehungen ersetzt:

∂r

∂τ= −vr (4.9)

∂ri

∂τ= −vi (4.10)

∂ri

∂τ=

∂

∂τ

(ri

r

)=

rivr − vi

r(4.11)

∂Un

∂τ=

∂Ui

∂τni + Ui

∂ni

∂τ

= Un + Un (4.12)

Angewandt auf die Gleichung von FW-H ergibt sich

4πp′(~x, t) =

∫

S

[%0(Un + Un)

r(1−Mr)2

]

τ=τ?

dS +

∫

S

[%0UnK

r2(1−Mr)3

]

τ=τ?

dS

+1

c

∫

S

[Lr

r(1−Mr)2

]

τ=τ?

dS +

∫

S

[Lr − LM

r2(1−Mr)2

]

τ=τ?

dS +1

c

∫

S

[LrK

r2(1−Mr)3

]

τ=τ?

dS

(4.13)

mit

K = rMr + cMr − cM2 (4.14)

LM = LiMi. (4.15)

Im Programm C2Noise werden die Quellterme auf den durchlassigen Oberflachen nach diesenFormeln berechnet. Um weitere Vereinfachungen vorzunehmen zu konnen, bewegen sich derBeobachter und der umstromte Korper mit der gleichen Geschwindigkeit. In diesem Fall istMr = 0, da die Machzahl fur jeden Quellpunkt in Richtung Beobachter konstant ist. DerTerm Un ist ebenfalls gleich null, da die Oberflache starr ist und die Ableitung des normiertenNomalenvektors somit null ist.Alle Geschwindigkeiten (Fluid und des Korpers) mussen zu allen Zeiten und an allen Ortenkleiner als die Schallgeschwindigkeiten sein, da es sonst zu Singularitaten bei der Berechnungder Quellterme kommt. Die Berechnung der retardierten Zeiten setzt ebenfalls voraus, daßzu allen Zeiten die Oberflache und der Beobachter sich langsamer als der Schall bewegen.


4.3.2 Die Bestimmung der retardierten Zeit

Die Bestimmung der retardierten Zeiten wird schon bei einfachen Bahnen der umstromtenKorper bei bewegtem Beobachter sehr aufwendig. Bei Unterschallgeschwindigkeit lautet dieBestimmungsgleichung fur die retardierte Zeit

c · (t− τ ?(~η)) =| ~x− ~η − ~xs(~η, τ ?(~η)) | (4.16)

mit ~η−~xs, dem momentanen Ortsvektor des Quellpunktes (3.70). ~η ist dabei der Ortsvektordes Quellpunktes zur Zeit τ = τ0.Eine exakte Gleichung zur Bestimmung der retardierten Zeiten laßt sich durch weitere Ver-einfachungen herleiten. Zur Zeit ist das Programm auf rein translatorische Bewegungen vonKorper und Beobachter in x1-Richtung, d.h. in Hauptstromungsrichtung, beschrankt.

v v

r

r

0

Quelle bei

t = τ0

Quelle beit = τ = τ + ∆τ0

∗

Beobachterx 1

v ∆τ.

Abbildung 14: Bestimmung der retardierten Zeit bei reiner Translation

Die Geschwindigkeiten konnen dabei unabhangig voneinander gewahlt werden und sind po-sitiv in die negative x1-Richtung. Unter diesen Bedingungen lassen sich die retardiertenZeiten auch analytisch bestimmen. Der Vektor ~η ist der Ortsvektor der Quellpunktes zurZeit τ = τ0, bei dem die Koordinatensysteme kongruent lagen.Dann gilt

τ =| ~x− ~η + v~ex1(τ − τ0) |

c. (4.17)

Bildet man aus Gleichung (4.17) eine quadratische Form in τ und ersetzt τ = τ0 +∆τ analogzur Abbildung (14), ergibt sich eine quadratische Formel in ∆τ :

∆τ 2 − 2∆τc2(t− τ0) + v(x1 − η1)

c2 − v2︸︷︷︸p

− r20 − c2(t− τ0)

c2 − v2︸︷︷︸q

= 0 (4.18)


mit der allgemeinen Losung

∆τ1/2 = p±√

p2 + q. (4.19)

Die retardierte Zeit ist dann mit folgender Gleichung bestimmbar:

τ ? = τ0 + p−√

p2 + q (4.20)

Gleichung (4.20) ist in der Funktion findtau zur Bestimmung der retardierten Zeiten imple-mentiert.

Es ist moglich, Gleichung (4.16) fur die iterative Berechnung der Quellzeit fur eine festeBeobachterzeit und einen festen Beobachterort umzuformen. Mit der Annahme, daß derBeobachter und der umstromte Korper sich geradlinig mit ~vB bzw ~v in beliebige Richtungenbewegen, gilt

c · (t− τ ?(~η)) =| ~x0 + ~vBt− ~η − ~v · τ ?(~η)) | . (4.21)

Wird die Gleichung iterativ mit dem Index k und dem Startwert τ ?0 ausgewertet, gilt fur

k ≥ 0 :

τ ?k+1 = t− | ~x0 + ~vBt− ~η − ~v · τ ?

k ) |c

(4.22)

Diese Formulierung hat einige Vorteile: Die Bewegungen von Beobachter und Quelle sindnicht auf die x1-Richtung beschrankt, die Berechnung der retardierten Zeiten benotigt je-doch mehr Rechenzeit, als die analytische Losung. Die Berechnung der retardierte Zeit kannbei Erreichen einer gewunschten Genauigkeit abgebrochen werden (sinnvoller Weise in derGroßenordnung der Zeitschrittweite der Stromungssimulation).Die iterative Berechnung der Quellzeitpunkte ist getestet worden, in der jetzt lauffahigenVariante von C2Noise wird das Verfahren jedoch nicht verwendet.

5.1 Konfiguration 45

5 Zweidimensionaler δ-Puls

Die Haupteffekte einer zweidimensionalen Schallquelle im Vergleich zu dreidimensionalenQuellen sind der erhohte Pegel und unendlich lange Signale bei kurzen Pulsen. Die zweidi-mensionalen Effekte werden auch bei großer endlicher Ausdehnung in x3-Richtung auftreten.Zur Bewertung dieses Effektes wird ein δ-Puls untersucht.

5.1 Konfiguration

Der Puls wird auf einem einzelnen Segment eines Kreises mit dem Durchmesser D = 0.1824merzeugt und mit verschieden großen Ausdehnungen in x3-Richtung der zweidimensionalenQuelle angenahert. Der Normalenvektor auf der Flache zeigt in negative x1-Richtung. DasKreissegment ist in Abbildung 15 dargestellt.

-100 -0.1 0.1

x1[m]

-0.1

0.10.1

x2[m]

nBeobachter

aktives Element i mit δ-Pulsn - Normalenvektor des aktiven Elementes

Abbildung 15: Element mit δ-Puls

Auf allen anderen Elementen ist der Druck fur alle Zeiten gleich null. Die Oberflache istundurchlassig. Der Testfall wird fur unbewegte Oberflache und Beobachter gerechnet. DerBeobachter befindet sich in 100m Entfernung zur Quelle bei den Koordinaten (x1, x2) =(−100, 0).Da die Quellterme des FW-H entweder eine Zeit- oder Ortsableitung enthalten, kann derδ-Puls nicht einfach auf der Oberflache vorgegeben werden. Der Puls wird durch einenDrucksprung der Form h(t; a) erzeugt. Im Intervall [−a, a] erfolgt der Anstieg von nullauf eins. Das Signal ist in Abbildung 16a dargestellt, der Drucksprung erfolgt dabei in-nerhalb von 100 Zeitschritten zum Zeitpunkt t0 = 0.5s. Die Pulsdauer ergibt sich damit zu∆TPuls = 0.01. Das komplette Signal ist definiert von 0.0001s bis 1.0s mit einer Zeitschritt-weite von ∆t = 0.0001s.

46 5.1 Konfiguration

h(t; a) =

0 fur t < −a12

(ta

+ 1π

sin(π ta) + 1

)fur −a ≤ t ≤ a

1 fur t > a

(5.1)

0.49 0.495 t0 0.505 0.51

t[s]

0

0.2

0.4

0.6

0.8

1

t0-a t

0t0+a

t[s]

0

0.2

0.4

0.6

0.8

1

h(t;a)

t0-a

t0+a

ddt__ h(t;a)

Abbildung 16: a) Signalverlauf und b) Ableitung des Signals

Der Quellterm auf einer festen Oberflache ist durch Gleichung (3.85) gegeben:

q(~x, t) =∂

∂t

{%0vn | gradf | δ(f)

}− ∂

∂xi

{Li | gradf | δ(f)

}(5.2)

mit Li = Pijnj = (p− p0)ni − τijnj.

Das Signal nach Gleichung (5.1) wird fur p eingesetzt. Da sich die Quelle nicht bewegt,fallt der Term mit der Zeitableitung heraus (vn ≡ 0 ∀t). Wird die Ortsableitung analog zuKapitel 4.3.1 umgeformt, ergibt sich die Losung als Zeitableitung des Drucksprungs. τijnj

wird dabei vernachlassigt.

4πp′(~x, t) =∂

∂t

∫

S

[(p− p0)rn

cr(1−Mr)

]

τ=τ?

dS +

∫

S

[(p− p0)rn

r2(1−Mr)

]

τ=τ?

dS

︸︷︷︸vernachlaessigt, da im Fernfeld sehr klein

(5.3)

5.1 Konfiguration 47

Der zweite Term mit dem Gleichanteil ist aufgrund der Abhangigkeit von 1r2 im Fernfeld

vernachlassigbar klein. Die Ableitung des Drucksprungs ist gerade ein δ-Puls mit der Breite2a.

δ(t; a) =d

dth(t; a)

= cos2(π

2

t

a

)fur − a < t < a. (5.4)

Der δ-Puls ist in Abbildung 16b dargestellt. Wird der Puls auf nur ein Element beschrankt,d.h. eine dreidimensionale Punktquelle definiert, so wird der Beobachter im Fernfeld denPuls in der Form unverandert mit verringerter Amplitude registrieren. Wird die Quelle inx3-Richtung fortgesetzt, werden 2D-Effekte auftreten, d.h. die abfallende Flanke des δ-Pulseswird sich verandern.

In den Abbildungen 17 und 18 sind die berechneten Signale des δ-Pulses beim Beobach-ter fur verschiedene Ausdehnungen in x3-Richtung dargestellt. Die Gesamtlange lz der 2D-Punktquelle uberdeckt den Bereich von ca. 4m bis 512m. Der Beobachter befindet sich in100m Entfernung. Der mit uber 500m langste Testfall approximiert ein zweidimensionaleQuelle ausreichend genau. Die einzelnen Testfalle sind in der folgenden Tabelle aufgelistet,lz ist dabei die Gesamtausdehnung der Quelle in x3-Richtung, nz die Anzahl und dz die Ele-mentbreite der Diskretisierung in x3-Richtung. In der Tabelle ist auch vermerkt, in welcherder Abbildungen die Testfalle dargestellt sind.

TPuls [s] Testfallname ] Elemente nz ∆z[m] lz[m] Abb. 17 Abb. 18 Abb. 190.01 n050.dz04 50 0.04 4.04 •0.01 n100.dz04 100 0.04 8.04 •0.01 n100.dz08 100 0.08 16.08 •0.01 n200.dz08 200 0.08 32.08 •0.01 n200.dz16 200 0.16 64.16 • •0.01 n200.dz32 200 0.32 128.32 • •0.01 n400.dz32 400 0.32 256.32 •0.01 n800.dz32 800 0.32 512.32 •0.002 n050.dz04 50 0.04 4.04 •0.002 n050.dz08 50 0.08 8.08 •0.002 n050.dz16 50 0.16 16.16 •0.002 n050.dz32 50 0.32 32.32 •0.002 n050.dz64 50 0.64 64.64 •

48 5.2 Auswertung

5.2 Auswertung

0.79 0.8 0.81 0.82t [s]

0

0.2

0.4

0.6

0.8

1

p’(t

)

n050.dz04n100.dz04n100.dz08n200.dz08n200.dz16n200.dz32

2D-δ-Puls mit TPuls

=0.01s

lz=32m

lz=64m

Abbildung 17: Zweidimensionaler δ-Puls fur verschiedene Langen

In Abbildung 17 ist erkennbar, wie der Pegel des Pulses mit großer werdender Ausdehnungin x3-Richtung anwachst. Eine Verdopplung der Lange lz verdoppelt nahezu auch den Pegel.Dies ist durch die Verdopplung der aktiven Flache erklarbar. Bei einer Pulsdauer von TPuls =0.01s wird der Maximalpegel bei einer Lange von lz = 64m erreicht.Wird die Dauer des Pulses auf ein Funftel verkurzt , TPuls = 0.002s, so wird der Maximal-pegel des Signals beim Beobachter schon bei geringerer Ausdehnung lz = 32m erreicht. DieTestfalle des kurzeren Pulses sind in Abbildung 19 dargestellt.

Da der lange δ-Puls mit 100 Zeitschritten aufgelost wird und in diesen Zeitraum von T =0.01s auf allen Elementen das Signal gleichphasig ist, addieren sich die Quellterme optimal,solange die aktiven Elemente zum gleichen Beobachterzeitintervall einen Beitrag liefern. Esgibt dann keine Terme, die fur eine Ausloschung des Signals sorgen.

Fur ein harmonisches Zeitsignal gilt dann:Je kleiner die Frequenz, desto weiter muß die Rechnung in x3-Richtung ausgedehnt werden,um den maximalen Pegel zu erreichen. Bei Variation der Ausdehnung des Korpers wirdbesonders bei tiefen Frequenzen der Unterschied im Pegel sehr groß sein.

5.2 Auswertung 49

0.8 0.9 1 1.1 1.2t [s]

0

0.2

0.4

0.6

0.8

1p’

(t)

n200.dz16n200.dz32n400.dz32n800.dz32


=0.01s

Abbildung 18: Zweidimensionaler δ-Puls fur verschiedene große Langen

Die Ursache dafur, daß der Pegel nicht immer weiter ansteigt, ist der zunehmende Winkelzwischen dem Normalenvektor und der Richtung zum Beobachter. Je weiter die Elementein x3-Richtung vom Ursprung entfernt sind, desto großer wird der Winkel zwischen ~r, demVektor vom Element zum Beobachter und ~n, dem Normalenvektor auf dem Element. Da-mit wird der effektive Wert der Quelle verringert, da im Quellterm aus Gleichung 5.3 dasSkalarprodukt von ~r und ~n, die Große ~rn steht. Der Anstieg des Pegels mit zunehmenderAusdehnung wird immer geringer. Dies ist gut bei den Testfallen n200.dz08 und n200.dz16

erkennbar. Bei weiterer Verlangerung der Quelle steigt der Pegel nicht weiter an.

Wird die Ausdehnung der Quelle weiter gesteigert, so bildet sich die Signalschleppe einerzweidimensionalen Schallquelle aus. In Abbildung 18 sind die Ausdehnung der Quelle von64m bis zu 512m dargestellt. Die Uberlagerung der δ-Pulse der weit außen liegenden Ele-mente ergibt den sanften Abfall des Signals. Der nur 0.01s lange δ-Puls erzeugt bei einerAusdehnung von ca. 500m der Quelle in x3-Richtung ein Signal von uber 0.4s. Das Ende desSignals lauft bei allen Langen immer sanft auf den Nullpegel aus.

50 5.3 Erforderliche Diskretisierung in x3-Richtung

0.794 0.796 0.798 0.8t [s]

0

0.2

0.4

0.6

0.8

1

p’(t

)

n050.dz04n050.dz08n050.dz16n050.dz32n050.dz64


=0.002s

lz=16m

lz=32m

lz=64m

Abbildung 19: Kurzer δ-Puls fur verschiedene Langen lz

5.3 Erforderliche Diskretisierung in x3-Richtung

Die zweidimensionalen Quelle, auf einem beliebigen Element in der x1 − x2−Ebene, wird inC2Noise durch identisches Fortsetzen der Quelle in der dritten Dimension approximiert. DieDiskretisierung ∆z der Elemente in x3-Richtung darf nicht beliebig grob werden, da ein Signalim Außenbereich eine andere Laufzeit zum Beobachter besitzt als beispielsweise ein Elementbei x3 = 0. Damit werden die Signale der Quellen im Außenbereich spater vom Beobachterregistriert, als Signale von Quellen, die naher an der x1−x2−Ebene liegen. Eine hinreichendfeine Diskretisierung ∆z in x3-Richtung soll gewahrleisten, daß die zweidimensionale Quellephysikalisch korrekt simuliert wird, indem die Elemente mit identischen x1−x2−Koordinatenden gleichen Quellwert besitzen, jedoch durch unterschiedliche x3 = j ·∆z mit j = 0, . . . , nzunterschiedliche retardierte Zeiten besitzen.

Fur eine Abschatzung der Diskretisierung ∆z werden fester Beobachter und Quelle voraus-gesetzt. Die Laufzeiten der Signale fur feste Beobachter- und Quellzeit sind dann nur vomAbstand r eines Elementes zum Beobachter abhangig.

Ein Element mit der Breite ∆x in x1-Richtung, mit identischem Quellwert in x3-Richtung,soll in seiner kompletten Ausdehnung ∆z die gleiche retardierte Zeit besitzen. Alle Punkteinnerhalb des Ringsegmentes mit der Breite c · ∆t besitzen die gleiche retardierte Zeit, dainnerhalb dieses Zeitraumes alle Punkte einen Beitrag zum selben diskreten Beobachterzeit-punkt liefern. Auf dem Kreis mit dem Radius r = c · (t− τ ?) um den Beobachter liegen alle

5.3 Erforderliche Diskretisierung in x3-Richtung 51

Punkte, die zu einem festen Zeitpunkt vom Beobachter registriert werden. Dieser Beobach-terzeitpunkt ist somit auch identisch mit dem Beobachterzeitpunkt fur alle Punkte innerhalbdes Ringsegmentes.

Der Fall ist in Abbildung 20 dargestellt. Es ist die x1 − x3−Ebene dargestellt.

R

R

r

r

x∆

z∆x3

min

z∆x3+

θ

Element i

x3

x1

c t.∆

Abbildung 20: Herleitung der maximalen Ausdehnung eines Elementes in x3-Richtung

R ist der Abstand des Beobachters zur Quelle in der x1 − x2−Ebene. r ist der Abstand imdreidimensionalen Raum.

Mit der Bedingung, daß die Breite des Elementes ∆x ≤ c · ∆t2

ist, ist garantiert, das solangeder Kreis mit dem Radius r = c · (t − τ ?) innerhalb des Elementes verlauft, alle Punktedes Elementes dem selben diskreten Beobachterzeitpunkt zugeordnet werden. Somit ist dieAnnahme des konstanten Quellwertes uber die gesamte Flache des Elementes mit der Lange∆z korrekt.

Die maximale Lange ∆zmax des Elementes i ist bestimmt durch die Schnittpunkte des Kreisesmit dem Radius r und den Grenzen des Elementes in der x1-Richtung. Die Ausdehnungdes Elementes i ist somit abhangig von ∆x, dem 2D-Abstand R und dem Abstand r desElementes zum Beobachter.

Aus Abbildung 20 sind zwei Beziehungen fur die Schnittpunkte von r mit dem Elementersichtlich.


r2 =(R +

∆x

2

)2

+ x23 und r2 =

(R− ∆x

2

)2

+(x3 + ∆z

)2. (5.5)

Das Gleichsetzen der beiden Gleichungen fuhrt auf eine quadratische Form in ∆z.

0 = ∆z2 + 2x3 ·∆z − 2R∆x (5.6)

Durch die Symmetrie des Problems unterscheiden sich die beiden moglichen Losungen nurdurch das Vorzeichen. Fur die maximale Ausdehnug des Elementes ergibt sich

∆z ≤√

x23 + 2R∆x− x3. (5.7)

Die Losung ist in Abbildung 21 geplotet.

0 50 100 150Ausdehnung in x

3-Richtung [m]

0.01

0.1

1

max

imal

e E

lem

entg

röß

e ∆z [m

]

x = 0.011m∆

R = 100m

Abbildung 21: Maximale Diskretisierung eines Elementes in x3-Richtung

Je langer der Korper in x3-Richtung ausgedehnt wird, desto kleiner mussen die Elementewerden. Ist der Beobachter weiter entfernt als in x3-Richtung gerechnet wird, kann die Dis-kretisierung ∆z großer als ∆x gewahlt werden. Wird die Rechnung weiter ausgedehnt, soist bei x3 = R die maximale Diskretisierung ∆z gleich ∆x. Bei weiterer Vergroßerung derAusdehnung muß ∆z kleiner als ∆x gewahlt werden.

Die Herleitung der Formel (5.6) erfolgte fur einen festen Beobachterzeitpunkt. Das Programmwertet jedoch immer Zeitintervalle ∆t aus und ordnet alle Signale der Quellen innerhalb c·∆teinem bestimmten Beobachterzeitpunkt zu. ∆zmax ist somit auch von der Zeitschrittweite

5.3 Erforderliche Diskretisierung in x3-Richtung 53

abhangig. Formel (5.6) gilt nur sofern c · ∆t2

großer als ∆x ist und sich der Quellwert desElementes sich zwischen zwei Quellzeitpunkten andert. Der Effekt wird also nur bei hoch-frequenten Anteilen der Losung auftreten.

Eine Betrachtung der maximalen Diskretisierung ∆z fur sehr große Ausdehnung der Quellein x3-Richtung ergibt, das ∆z, analog zu ∆x, nicht kleiner werden muß als ∆z = c · ∆t

2. Wird

∆z so groß gewahlt, ist immer gewahrleistet, das das komplette Element in jeder Richtunginnerhalb eines diskreten Quellzeitintervalles liegt, und damit korrekt einer retardierten Zeitzugeordnet wird.

Anhand des δ-Pulses soll nun noch gezeigt werden, wie sich eine zu grobe Diskretisierung∆z auf die Simulationsergebnisse auswirkt.In Abbildung 22 sind die Signale an der Beobachterposition fur einen δ-Puls mit der BreiteTPuls = 0.002s dargestellt. Der Beobachter befindet sich wieder in 100m Entfernung ander Position (x1, x2) = (0, 100) und die beiden Testfalle sind ca. 128m lang, jedoch mitunterschiedlichen Auflosungen ∆z gerechnet. Die Parameter der Akustikrechnung sind infolgender Tabelle aufgelistet.

Testfallname ] Elemente nz ∆z[m] lz[m] Abbildung 22a Abbildung 22bn400.dz16 400 0.16 128.16 • •n100.dz64 100 0.64 128.64 • •


0.78 0.8 0.82 0.84t [s]

0

0.2

0.4

0.6

0.8

1p’(t) n400.dz16n100.dz64

0.83 0.831 0.832 0.833 0.834t [s]

0.07

0.08

0.09

0.1

p’(t)a) b)

Abbildung 22: Signal eines zweidimensionalen δ-Pulses mit unterschiedlicher Diskretisierung ∆zgerechnet: a) kompletter Puls, b) Ausschnittvergroßerung

Wird mit zu grober Auflosung ∆z gerechnet, kommt es zu Schwingungen in der abfallendenFlanke des Signals. Die Schwingungen durch die zu grobe Auflosung sind in Abbildung22b vergroßert dargestellt. Eine weitere Steigerung der Diskretisierung fuhrt zu hoherenAmplituden der Schwingungen, andert jedoch nichts an der Frequenz.Bei zu grober Diskretisierung wird der Quellterm eines Elementes durch die zu große Flacheuberschatzt, und das beobachtete Signal zu diesem Zeitpunkt zu groß sein. Das Signal furden Zeitschritt ∆t danach wird um den gleichen Wert zu niedrig vorhergesagt.

Als Abschatzung des Dikretisierungsfehlers werden die Zeitsignale einer Simulation mit gro-ber Auflosung (n100.dz64) von einer mit feiner Auflosung (n400.dz16) subtrahiert. DasDifferenzsignal wird auf seine Frequenzanteile untersucht. Das Ergebnis ∆p′(f) ist in Abbil-dung 23 dargestellt.Die Frequenz mit dem großten Pegel liegt bei 1175Hz. Die ersten drei Oberwellen dieserGrundfrequenz sind ebenfalls mit hohen Pegeln vertreten. Dies korrospondiert mit dem peri-odischen An- und Abschwellen der Amplitude der hochfrequenten Schwingung im Zeitsignaldes Testfalls n100.dz64.

5.4 Zusammenfassung 55

1000 2000 3000 4000 5000f [Hz]

0.5

1

∆p’ (f) Fouriertransformation des Differenzsignals der Testfälle n100.dz64 und n400.dz16

Abbildung 23: Frequenzdarstellung der Differenz der beiden Testfalle

5.4 Zusammenfassung

Die Implementation des Algorithmus erlaubt eine sehr gute Approximation einer zweidimen-sionalen Schallquelle. Die Wahl der Ausdehnung des Korpers in der dritten Dimension inder Akustik, hat entscheidenden Einfluß auf den berechneten Pegel. Der Fehler des vorher-gesagten Pegels gegenuber der unendlich langen zweidimensionalen Quelle ist bei gegebenerLange auch frequenzabhangig. Je hoher die Frequenz, desto geringer die Abweichung.Fur die Simulation des Schalles auf Basis von zweidimensionalen Stromungssimulationen vonrealen Korpern, die endlich in x3-Richtung sind, gilt dann die folgende Beschrankung:Der berechnete Gesamtpegel wird nur durch die Wahl der effektiven akustischen Lange lzbestimmt. Eine Aussage uber die Absolutpegel ist auf Basis des vorliegenden implementiertenAlgorithmus nicht moglich.

Die Diskretisierung ∆z darf nicht zu grob werden, da sonst die Ergebnisse im hohen Fre-quenzbereich stark fehlerbehaftet sind. Eine Formel zur Abschatzung ist durch Gleichung(5.6) gegeben.Die maximale Diskretisierung ∆z ist von der Entfernung des Beobachters zur Quelle abhangig,und sinkt mit zunehmender Ausdehnung der Quelle in x3-Richtung. Ist der Beobachter wei-ter von der Quelle entfernt als die Quelle in x3-Richtung ausgedehnt ist, so kann ∆z einVielfaches der Diskretisierung in der x1 − x2-Ebene betragen.Allgemein sollte die Diskretisierung ∆xi kleiner sein als c· ∆t

2. In der vorliegenden zweidimen-

sionalen Implementierung kann ∆z auf Basis der Gleichung (5.6) grober gewahlt werden, dadiese das obige Kriterium bei sehr großer Ausdehnung der Quelle automatisch erfullt.

56 6.1 Stromungssimulation

6 Zylinderumstromung

Die Simulation der Umstromung eines Kreiszylinders als erstem Testfall erfolgt aufgrundvorhandener Akustikmeßdaten von King III u.a. [6] an Zylindern mit verschiedenen Langen-verhaltnissen. Die zur Verfugung stehenden experimentellen Messdaten von King sind sehrausfuhrlich. Alle Parameter fur die Stromungssimulation sind dieser Arbeit entnommen.

Des weiteren ist eine Veroffentlichung von Cox u.a. [3] vorhanden, die sich mit Stromungssi-mulationen und anschließender Akustikberechnung auf Basis des FW-H eines Zylinders beiverschiedenen Re-Zahlen beschaftigt. Cox zeigt dabei den Einfluß verschiedener Turbulenz-modelle auf die Ergebnisse der Simulationen und vergleicht mit experimentell ermitteltenDaten.

Die Reynoldszahlen der Testfalle betragen Re = 86400 bei King III und Re = 89000 (expe-rimentelle Daten) bzw. Re = 90000 (Simulation) bei Cox.

6.1 Stromungssimulation

6.1.1 Konfiguration

Die zweidimensinale Stromungssimulation wird mit ELAN2 mit den folgenden Parameterngerechnet:

Der Durchmesser des Zylinders betragt D = 0.02m bei einer Anstromgeschwindigkeit vonu∞ = 64m

s. Dies entspricht einer Machzahl von M = 0.188. Die kinematische Viskositat

betragt ν = 1.52 · 10−5 m2

sund die Dichte % = 1.2 kg

m3 .

Damit ergibt sich die mit dem Durchmesser gebildete Reynoldszahl Re = 86400.

Die Stromungssimulation wird dimensionsbehaftet, vollturbulent und da die Machzahl kleinist, inkompressibel gerechnet. Der Turbulenzgrad der Grundstromung, vorgegeben durch dasExperiment, ist auf Tu = 0.5% eingestellt.

In zeitlicher Richtung wird eine Zeitschrittweite ∆t = 0.0001s bei 9990 Zeitschritten verwen-det. Die Gesamtsimulationszeit betragt ca. TA = 1s. Die Zeitdiskretisierung ist von zweiterOrdnung genau.

Turbulente Effekte werden durch das LLR-k − ω−Modell, ein modifiziertes k − ω-Zwei-gleichungsmodell, erfaßt. Die Konvektionsterme werden mit Hilfe des QUICK-Schemas dis-kretisiert.

6.1.2 Rechengitter

Das Gitter besteht aus 16068 Volumenelementen in drei Blocken. Die Volumenelementesind zur Zylinderoberflache stark verdichtet. Die Gesamtausdehnung betragt in der Hohe 30Zylinderdurchmesser und 48 Zylinderdurchmesser in Hauptstromungsrichtung. Der Zylinderist 8 Zylinderdurchmesser vom Einstromrand entfernt.

Abbildung 24 zeigt das Rechengitter in unmittelbarer Umgebung des Kreiszylinders, mit derAnstromung von links.

6.1 Stromungssimulation 57

Abbildung 24: Verwendetes Rechengitter um den Kreiszylinder

Das Gitter erlaubt das Rechnen mit LOW-REYNOLDS-Randbedingungen, d.h. die Grenz-schicht wird komplett bis in den linearen Bereich aufgelost. Die normierte Hohe y+ = yuτ

ν

der wandnachsten Zellen liegt im Bereich von 0.2 ≤ y+ ≤ 2.0.

6.1.3 Randbedingungen

Auf der Zylinderoberflache wird Wandhaftung angenommen. Das Fluid bewegt sich an derWand mit der gleichen Geschwindigkeit wie die Korperoberflache.Auf dem Einstromrand sind die Geschwindigkeiten und turbulenten Großen vorgegeben. DieGeschwindigkeit wird als konstant angenommen.Am Ausfluß wird eine konvektive Randbedingung verwendet, oben und unten wird Symme-trie gefordert.


6.1.4 Experimentelle Vergleichsdaten

Experimentelle Untersuchungen der Umstromung von Kreiszylindern ergeben im Bereich derReynoldszahlen unter Re = 1.0 · 105 eine gleichmaßige Wirbelablosung mit Strouhalzahlenum St = 0.195± 0.05.

2*104

105

106

107

ReD

0.1

0.2

0.3

0.4

0.5

St

exp. Strouhal-Zahlen eines Kreiszylinders (Schewe, 1983)

Abbildung 25: Strouhalzahl uber Re eines Kreiszylinders, experimentelle Werte nach Schewe [13]

In Abbildung 25 sind in Experimenten bestimmte Strouhalzahlen uber einem weiten Reynolds-zahlbereich aufgetragen.Die Umstromung des Zylinders bei Re = 86000 liegt im unterkritischen Bereich. Experimen-telle Werte fur den Widerstandsbeiwert cd nach Cox [3] liegen im Reynoldszahlbereich von50000 ≤ Re ≤ 1.0 · 105 bei 0.9 ≤ cd ≤ 1.4.Die kritische Reynoldszahl des Kreiszylinders liegt bei etwas uber Re = 1.0 ·105. Der Wider-standsbeiwert sinkt im Bereich von 1.0 ·105 ≤ Re ≤ 4.0 ·105 rapide auf rund cd = 0.3, da dieUmstromung turbulent wird. Der Auftriebsbeiwert cl des Zylinders liegt im Mittel bei null.


6.1.5 Auswertung

In Abbildung 26 ist der instationare Verlauf des Auftriebs- und Widerstandsbeiwertes derStromungssimulation des Kreiszylinders bei Re = 86400 dargestellt. Zum Vergleich sind inAbbildung 27 die Ergebnisse von Cox [3] fur verschiedene Turbulenzmodelle bei Re = 90000dargestellt.

0.3 0.302 0.304 0.306 0.308t [s]

-2

-1

0

1

2c

l

0.3 0.302 0.304 0.306 0.308t [s]

0

0.5

1

1.5

2c

d

Abbildung 26: Instationarer Verlauf des Auftriebsbeiwertes cl und Widerstandsbeiwertes cd einesKreiszylinders bei Re=86400

Abbildung 27: Instationarer Verlauf des Auftriebsbeiwertes cl und Widerstandsbeiwertes cd einesKreiszylinders bei Re = 90000, vollturbulente Stromungssimulation nach Cox [3]

Die Schwankungen von cl sind sehr gleichmaßig um den Mittelwert von cl = 0 herum, derZylinder erzeugt keinen Auftrieb. Die Amplituden liegen im Bereich von |cmax

l | ≈ 1.4. DieVergleichswerte von Cox liegen durchweg darunter, eine nicht dargestellte laminare Rechnungergibt auch bei Cox Amplituden von |cmax

l | = 1.4. Der Widerstandsbeiwert cd schwankt umden Mittelwert von cd = 1.255 zwischen |cmin

d | = 1.0 und |cmaxd | = 1.55. Die turbulenten

Rechnungen von Cox liegen durchweg darunter, die laminare Rechnung liegt ebenfalls imMittel uber eins.Die zeitlichen Verlaufe weisen bei Cox eine starke Abhangigkeit vom verwendeten Turbu-lenzmodell auf. Das S-A-Modell beispielsweise zeigt kaum instationare Effekte.


Die vorliegende Simulation ist vollturbulent gerechnet. Die großen Differenzen zu den tur-bulenten Simulationen von Cox sind wahrscheinlich auf den sehr kleinen Turbulenzgrad derAnstromung Tu = 0.5% zuruckzufuhren. Die Umstromung des Zylinders ist laminar undwird erst nach der Ablosung turbulent.Cox gibt nicht den verwendeten Turbulenzgrad und andere Randbedingungen an. Die guteUbereinstimmung mit den laminaren Rechnungen von Cox bestatigen jedoch die Vermutun-gen.

Die Umstromung des Zylinders bei Re = 86400 erfolgt laminar und ist durch abwechselndesregelmassiges Ablosen von Wirbeln an Ober- und Unterseite gekennzeichnet. Der Punkt andem die Stromung vom Korper ablost, schwankt sehr stark in der Zeit. Der Umschlag vonlaminar nach turbulent erfolgt erst kurz nach der Ablosung. Analog zur Ablosung schwanktauch die Position des Staupunktes mit der Zeit. Es bildet sich eine Wirbelstraße hinter demZylinder.Die Strouhalzahl, gebildet aus dem Zeitverlauf des Auftriebsbeiwertes cl betragt rund St =0.26. Die Strouhalzahl

St = f · l

u∞(6.1)

wird mit dem Durchmesser D = 0.02m, der Anstromgeschwindigkeit u∞ = 64ms

und derFrequenz von ca. f = 830Hz gebildet.Die Strouhalzahl von St = 0.26 liegt weit entfernt von den experimentellen Daten vonStexp ≈ 0.195. Die Vergleichsrechnungen von Cox bei Re = 90000 ergeben Strouhalzahlenim Bereich von 0.21 ≤ StCox ≤ 0.29. Das Turbulenzmodell hat entscheidenden Einfluß aufdie vorhergesagte Strouhalzahl, bei der laminaren Rechnung von Cox ergibt sich St = 0.235.

Verbesserungen der Ergebnisse durch Variation der Randbedingungen sind nicht festzustel-len. Eine Vergleichsrechnug mit gesetzter Transition im Bereich der Umstromung mit be-schleunigter Grenzschicht , d.h. vom Staupunkt bis zu Punkt xTransition = −0.08, verandertdie Strouhalzahl nicht, da die Umstromung bei der vollturbulenten Rechnung ebenso laminarist.

URANS-Rechnungen sind fur die Simulation der Umstromung von stumpfen Korpern naheder kritischen Reynoldszahl nicht gut geeignet. Die starken Schwankungen der Transitionund Ablosung konnen nur unzureichend von URANS-Rechnungen physikalisch korrekt wie-dergegeben werden. Bessere Ergebnisse sind nur mit LES (large eddy simulation) bzw. DES(detached eddy simulation) moglich, vergleiche Schmidt [14].

6.2 Akustiksimulation 61

6.2 Akustiksimulation

6.2.1 Experimentelle Vergleichsdaten

Die Akustikmessungen von King III u.a. [6] erfolgten an zwei Punkten, senkrecht zur An-stromung auf gegenuberliegenden Seiten des Zylinders. Die Messmikrofone sind bei King imAbstand von 1.4m positioniert.King untersucht dabei den Einfluß verschiedener Langenverhaltnisse L/D des Kreiszylindersauf den Schalldruckpegel.

Abbildung 28: Messung der Spektren eines Kreiszylinders D = 0.02m, uo = 64ms , senkrecht zur

Anstromung in 1.4m Abstand, fur verschiedene Langenverhaltnisse L/D nach King III [6]

In den Abbildungen 28 und 29 sind die Messungen fur die Langenverhaltnisse im Bereichvon 15 ≤ L/D ≤ 35 dargestellt. Die Bezeichnung des Testfalls A1 von King definiert einenZylinder mit dem Durchmesser D = 0.02m.King stellt Schmalbandspektren von 200Hz bis 1000Hz dar. Mit zunehmender Lange desZylinders steigen die Wirbelablosefrequenz und der Pegel dieser. Bei L/D = 15 betragt diePeakfrequenz fPeak = 586Hz entsprechend einer Strouhalzahl St = 0.183 und dem Peakpegelvon SPLPeak = 72dB. Bei L/D = 20 steigt die Frequenz auf fPeak = 602Hz bei einem Pegelvon SPLPeak = 85.5dB. Die gemessenen Signale sind sehr breitbandig, dreidimensionaleEffekte sind stark ausgepragt.

Ab einem Langenverhaltnis von L/D ≥ 25 ist keine Steigerung des Maximalpegels derPeakfrequenz mehr ersichtlich. Der Maximalpegel betragt rund SPLPeak = 90.5dB, bei einerFrequenz von fPeak = 616Hz. Dies entspricht einer Strouhalzahl von St = 0.193. Bei demLangenverhaltnis L/D = 25 ist der Frequenzpeak sehr schmal, bei weiterer Verlangerung desZylinders nehmen die dreidimensionalen Effekte wieder zu, das Spektrum wird breiter. DerSchalldruckpegel SPL auf Basis des gesamten Spektrums steigt im Gegensatz zum PeaklevelSPLPeak jedoch bei Vergroßerung der Lange auf L/D = 35 immer weiter an.Abschließend sind die Ergebnisse der Messungen von King in folgender Tabelle dargestellt.Der Schalldruckpegel SPL in dB basiert auf dem RMS-Werte, gebildet im Frequenzbereich

62 6.2 Akustiksimulation

Abbildung 29: Messung der Spektren eines Kreiszylinders D = 0.02m, uo = 64ms , senkrecht zur

Anstromung in 1.4m Abstand, fur verschiedene Langenverhaltnisse L/D nach King III [6]

von 100Hz bis 1000Hz. Er ist deutlich hoher als der Pegel der Peakfrequenz fPeak, da dieSignale sehr breitbandig sind.

L/D u∞[ms] SPL[dB] fPeak[Hz] St

3 64 81.6 - -9 64 84.0 540 0.16915 64 87.1 586 0.18317.5 64 90.1 598 0.18720 64 92.5 602 0.18825 64 98.0 616 0.19330 64 100.5 608 0.19035 64 102.3 608 0.190

Ausgepragte zweidimensionale Signale sind bei King mit Ausnahme des LangenverhaltnissesL/D = 25 nicht zu erkennen. Die dreidimensionalen Effekte dominieren die Akustikmessun-gen. Das Langenverhaltnisses von L/D = 25 dient als Vergleichsmesswerte fur die eigenenSimulationen.

Die Messwerte von Revell aus der Veroffentlichung von Cox u.a. [3] beschranken sich aufein Spektrum, gemessen 90 Grad zur Anstromung bei Re = 89000, Anstromgeschwindigkeitu∞ = 68m

sund einem Durchmesser D = 0.02m. Die Entfernung von Mikrofon zum Zylinder

betragt 2.56m. Die Strouhalzahlen 0.04 ≤ St ≤ 0.62 des Spektrums in Abbildung 30 ent-sprechen Frequenzen von 140Hz ≤ f ≤ 2100Hz. Die Peakfrequenz liegt bei St = 0.19. Diesentspricht einer Frequenz von fPeak = 646Hz.


In Abbildung 30 sind die Ergebnisse der Simulationen von Cox dargestellt.Diese werden zumVergleich der eigenen FW-H-Berechnungen genutzt.

Abbildung 30: Simulierte Spektren eines Kreiszylinders D = 0.02m, Re = 90000 senkrecht zurAnstromung in 2.56m Abstand, nach Cox [3] in Vergleich mit Experimenten von Revell bei Re =89000

Wahrend die Peakfrequenzen der beiden Messungen von Cox und King gut ubereinsimmen,unterscheiden sie sich signifikant im Pegel der Peakfrequenz. Bei Cox liegt der Pegel derPeakfrequenz bei SPLCox

Peak = 99dB, bei King liegt das Maximum bei SPLKingPeak = 90.5dB.

Der Unterschied von knapp 10dB im Pegel ist gravierend. Die Messung bei King erfolgte ingeringerer Entfernung als bei Cox, die Messung liefert jedoch einen geringeren Pegel. DerPegel sollte jedoch mit zunehmenden Abstand geringer werden.

6.2.2 Konfiguration

In Abbildung 31 sind die Kontrollflachen fi dargestellt. Auf diesen Flachen sind die in-stationaren Druckdaten und Geschwindigkeitskompomenten aus der Stromungssimulationextrahiert. Die Oberflache f1 ist identisch mit der Oberflache des Zylinders, hier sind nurdie Druckdaten vorhanden. Die anderen Kontrollflachen sind durchlassig. Die Flachen sindalle kreisformig und um den Zylindermittelpunkt angeordnet.Die Flachen f2 und f3 sind in unmittelbarer Nahe des Zylinders mit dem 1.5fachen bzw.doppelten Zylinderdurchmesser D positioniert. Flache f4 hat den Durchmesser d = 4D undf5 einen Durchmesser von d = 5D. Mit f4 und f5 soll untersucht werden ob sehr weit vonder Korperoberflache entfernte Kontrollflachen stark abweichende Ergebnisse produzieren.

Die Parameter fur die Berechnungen mit C2Noise sind durch die Stromungssimulation vor-gegeben: Die Kontrollflache und der Beobachter bewegen sich im Bezugssystem mit v = 64m

s

nach links. Der Schalldruckverlauf wird an den Messpositionen von King in rmess = 1.4m,90Grad zur Anstromung berechnet. Zusatzlich wird auf einem Kreis mit dem Radius r =


-0.05 0 0.05 0.1x [m]

-0.06

-0.04

-0.02

0

0.02

0.04

0.06

y [m

]

f1: undurchlässige Zylinderoberfläche

f2: durchl. Kontrollfläche, d = 1.5D

f3: durchl. Kontrollfläche, d = 2D



Abbildung 31: Kontrollflachen fi, auf denen die Druck- und Geschwindigkeitsdaten aus derStromungssimulation extrahiert werden

1.4m um den Zylinder herum an 100 gleichmaßig verteilten Punkten der Schalldruckver-lauf berechnet, um die Richtcharakteristik der Schallabstrahlung des Zylinders darstellen zukonnen.

Die wichtigsten Parameter sind die Diskretisierung in z-Richtung und die effektive akusti-sche Lange, d.h. wie weit die Rechnung in z-Richtung ausgedehnt wird. Uber diese Parameterwird der Maximalpegel bestimmt, der simuliert wird. Die effektive Gesamtlange lz fur den2D-FW-H muß kleiner sein als die reale Ausdehnung des Zylinders von l = 0.5m, da dieUmstromung des Zylinders nicht auf der gesamten Lange gleichmassig erfolgt. Zwei Punk-te auf der Zylinderoberflache, die sich nur in der z-Koordinate unterscheiden, werden mitzunehmenden Abstand unkorrelierte Werte des Druckes auf der Oberflache aufweisen. Beigroßem Abstand konnen die Signale gegenphasig sein. Untersuchungen von Schmidt [14] aufBasis von DES-Simulation ergeben bei Reynodszahlen um Re = 1.4 · 105 schon im Abstandvon z1 − z2 = 3D stark unkorrelierte Signale.

Fur den Zylinder sind die optimale Parameter fur C2Noise die Anzahl der Elemente inz-Richtung nz = 10 bei einer Auflosung von ∆z = 0.001m. Es ergibt sich die effektiveLange lz = 0.021m. Bei vorliegender Wahl der effektiven akustischen Lange betragt derSchalldruckpegel senkrecht zur Anstromung ca. 93dB. Damit liegt der Pegel der Simulationzwischen denen der beiden Messungen.

6.2.3 Auswertung der FW-H-Simulation

Alle folgenden Auswertungen sind fur Beobachter im Abstand von 1.4m vom Zylindermit-telpunkt. Die Gradangaben in den Abbildungen beziehen sich auf den Winkel zwischen derHauptstromungsrichtung und der Richtung zum Beobachter. Die Angaben θ = 90o undθ = 270o sind senkrecht zur Anstromung und stimmen mit den Meßpositionen von Kinguberein. θ = 180o ist entgegen der Anstromung des Zylinders.


6.2.4 FW-H auf der festen Korperoberflache f1 ausgewertet

In Abbildung 32 ist ein Auschnitt des Zeitverlaufes der Schalldruckschwankungen p′ derBeobachter im Winkel θ = 90o und θ = 270o dargestellt.

0.309 0.31 0.311 0.312 0.313 0.314 0.315 0.316t [s]

-2

-1

0

1

2

p’(t

) [P

a]90

o

270o

- Beobachter in 1.4m Entferung vom Zylinder

- Testfall n10.dz001, feste Oberfläche

Abbildung 32: Zeitverlauf des Schalldruckes beim Beobachter im Abstand r = 1.4m

Der Signalverlauf an den gegenuberliegenden Beobachterpositionen ist gleich, jedoch um180o phasenverschoben. Dies korrospondiert mit dem wechselseitigen Ablosen der Wirbelvom Zylinder.Die Periodendauer des Signals beim Beobachter ist identisch mit der Schwankung des Auf-triebsbeiwertes (siehe Abbildung 26a). Die Periodendauer betragt rund T = 0.001196s.Der Signalverlauf des Schalldruckes p′ beim Beobachter (Abbildung 32)ist jedoch nicht sogleichmaßig wie der des Auftriebsbeiwertes. In den an- bzw. absteigenden Flanken des Si-gnals ist der Einfluß von hoheren Harmonischen, d.h. Oberschwingungen mit Frequenzen desMehrfachen der Grundfrequenz zu erkennen. Die Grundfrequenz wird als erste Harmonischebezeichnet, Tone mit der doppelten Frequenz als zweite Harmonische usw.Der Einfluß der hoheren Harmonischen außert sich in einem ’s’-Schlag im Signalverlauf.

Die Frequenzdarstellung des 90o-Signals, berechnet mittels einer diskreten Fouriertransfor-mation mit dem Programm XmGrace, ist in Abbildung 33 dargestellt. Die Ordinate gibt denSchalldruckpegel SPL in dB relativ zum Bezugspegel p0 = 2 · 10−5Pa an. Die Fouriertrans-formierte ist mit einem Hanningfenster berechnet, um Fehler durch die Periodisierung desSignals klein zu halten.In Abbildung 33 sind ebenfalls die experimentellen Daten von King III fur ein Langen-verhaltnis L/D = 25 aus Abbildung 29 und von Revell aus der Veroffentlichung von Cox ausAbbildung 30 dargestellt.Die Grundfrequenz bei fPeak = 836Hz mit einem Pegel von SPLPeak = 96dB. Die zweiteHarmonische ist mit einem Pegel von SPL2.H = 74dB und die dritte Harmonische miteinem Pegel von SPL3.H = 77dB vertreten. Sie sind damit ca. 22dB bzw. 19dB leiser als derGrundton.


200 400 800 2k1k 3kf [Hz]

0

10

20

30

40

50

60

70

80

90

100

110

SP

L [d

B, r

el. 2

*10-5

Pa]

LLR-k-ω/FW-H-Simulation: feste Zylinderoberfläche,

nz010.dz001, Re=86400, Abstand r=1.4m, 90o

Messung King III, Re=86400, Abstand r=1.4m, 90o

Messung Revell, Re=89000, Abstand r=2.56m, 90o

Abbildung 33: Spektrum beim Beobachter auf Basis der 2D-Stromungssimulation und FW-Hsenkrecht zur Anstromung des Zylinders in 1.4m Entfernung im Vergleich mit den Messungen vonKing III [6] und Revell [3]

Die Peaks sind sehr schmal, was auf die zweidimensionale CFD/FW-H zuruckzufuhren ist.Alle dreidimensionalen Effekte in der Stromung und Akustik werden komplett vernachlassigt.Diese sind fur die breitbandigen Signale der Messwerte verantwortlich.

Im Vergleich mit den Messwerten von King und Cox liegt die Frequenz des Grundtons aufBasis der FW-H-Simulation hoher. Dies korrespondiert mit den Unterschieden der Wirbe-lablosefrequenz von Messung und Stromungssimulation. Die Messwerte von Revell aus derVeroffentlichung von Cox in Abbildung 33 zeigen auch die Maxima der hoheren Harmoni-schen. Die zweite und dritte Harmonische der Messungen liegen ca. 27dB unter dem Pegeldes Grundtons. Dieses Verhaltnis der Pegel von Grund- zu Obertonen wird auch von derSimulation gut wiedergegeben.

Die Messwerte von King bieten in diesem Punkt keine Vergleichsmoglichkeiten, da die dar-gestellten Spektren zu schmalbandig sind (f ≤ 1000Hz).

Die Simulationen von Cox, auf Basis von Stromungssimulation und FW-H analog zu demAnsatz in dieser Arbeit, berechnen die Frequenzen mit ahnlich grossem Fehler. Die Peakfre-quenz schwankt nach Cox stark in Pegel und Frequenz, je nach genutzten Turbulenzmodellbei der Stromungssimulation (siehe Abbildung 30).

Die Qualitat der Stromungssimulation hat entscheidenden Einfluß auf die Akustiksimulationmit dem FW-H-Algorithmus.

Die Simulation der Schallausbreitung ermoglicht die Berechnung des Schalls an Stellen, dieim akustischen Windkanal nicht oder nur sehr schwierig durchfuhrbar sind. An 100 Stellenauf einem Kreis im Abstand von 1.4m um den Zylinder wird der Schalldruck berechnet, um


die Richtcharakteristik der Schallabstrahlung des Zylinders zu erhalten. Der ’Directivity-Plot’ ist in Abbildung 34 dargestellt.

100 95 90 85 85 90 95 100

SPL [dB, rel. 2*10-5

Pa]100

95

90

85

85

90

95

100

Zylinder: - auf fester Oberfläche

- TA = 1.0s

ausgewertet

Abbildung 34: Richtcharakteristik des Zylinders im Abstand von r = 1.4m bei Re = 86400

Der Schalldruckpegel basiert auf dem RMS-Wert, gebildet aus dem ca. TA = 1s langenZeitsignal am Beobachterort. Die Formeln zur Berechnug sind durch Gleichung 3.2 bis 3.3gegeben.Der Schalldruckpegel liegt im Bereich von 82.25dB ≤ SPL ≤ 93.75dB. Die Richtcharakte-ristik ist ahnlich dem eines Dipols. Der abgestrahlte Schall ist in Bewegungsrichtung kleinerals senkrecht zur Anstromung, gut sichtbar durch die Einschnurung im Directivity-Plot. DerPegel hinter dem Zylinder ist kleiner als in Richtung der Bewegung.Der maximale Schalldruckpegel liegt nicht in 90o bzw. 270o Richtung, sondern bei ca. 97o/263o

in Richtung der Anstromung. Analog zum bewegten Dipol steigt der Pegel in Bewegungs-richtung. Der Winkel ϑmax mit dem maximalen Pegel ist beim bewegten Dipol durch dieBeziehung sin ϑmax = M gegeben. Er ist senkrecht zur Bewegungsrichtung definiert. Fur denbewegten Dipol betragt der Winkel bei einer Machzahl von M = 0.188 ϑmax = 10.84o. BeimZylinder ist der Winkel kleiner. Das Signal enthalt im Gegensatz zum idealen Dipol nichtnur eine Frequenz. Die Pegel der Harmonischen hat ebenso Einfluß auf den Gesamtpegel.


Die Zeitsignale an den verschiedenen Positionen geben einen Anhaltspunkt, weshalb derPegel in Richtung der Anstromung sinkt. In Abbildung 35 sind die Zeitsignale senkrecht zurAnstromung und im Winkel θ = 111.6o dargestellt.

0.309 0.31 0.311 0.312 0.313 0.314 0.315 0.316t [s]

-2

-1

0

1

2

p’(t

) [P

a]

0.309 0.31 0.311 0.312 0.313 0.314 0.315 0.316t[s]

-2

-1

0

1

2

p’(t

) [P

a]

90o

- Testfall n10.dz001, feste Oberfläche, Beobachter in 1.4m Entfernung

111.6o


Der ’s’-Schlag im Zeitsignal wird mit zunehmenden Winkel starker. Die zweite Harmonischewird im Verhaltnis zum Grundton lauter. Besser ist dies in den Spektren der entsprechendenZeitsignale zu erkennen. Sie sind in Abbildung 36 dargestellt.

200 400 800 2k1k 3kf [Hz]

0

20

40

60

80

100

SP

L [d

B, r

el. 2

*10-5

Pa]

200 400 800 1k 2k 3k400 800 2k1k 3kf [Hz]

0

20

40

60

80

100

SP

L [d

B, r

el. 2

*10-5

Pa]

90o

- Testfall n10.dz001, TA = 1.0s, feste Oberfläche, Beobachter in 1.4m Entfernung

111.6o

Abbildung 36: Spektrum des Schalldruckes beim Beobachter im Abstand r = 1.4m

Der Pegel des Grundtons und der dritten Harmonischen bleiben fast konstant. Bei der Ver-grosserung des Beobachtungswinkels von θ = 90o auf θ = 111.6o steigt die zweite Har-monische um 6dB, verdoppelt also die Lautstarke. Wird der Beobachtungswinkel weitervergrossert, steigt der Einfluß der hoheren Harmonischen weiter an. Dies korrespondiert mitder zunehmenden Starke des ’s’-Schlags im Signal.


0.309 0.31 0.311 0.312 0.313 0.314 0.315 0.316t [s]

-2

-1

0

1

2

p’(t

) [P

a]

0.309 0.31 0.311 0.312 0.313 0.314 0.315 0.316t[s]

-2

-1

0

1

2

p’(t

) [P

a]


136.8o


158.4o


In Abbildung 37 sind Auschnitte aus dem Schalldruckverlauf in den Richtungen θ = 136.8o

und θ = 158.4o dargestellt. Die entsprechenden Spektren sind in Abbildung 38 dargestellt.

200 400 800 2k1k 3kf [Hz]

0

20

40

60

80

100

SP

L [d

B, r

el. 2

*10-5

Pa]

200 400 800 1k 2k 3k400 800 2k1k 3kf [Hz]

0

20

40

60

80

100

SP

L [d

B, r

el. 2

*10-5

Pa]

136.8o


158.4o


Der Pegel der zweiten Harmonischen steigt mit zunehmenden Winkel um je 4dB bzw. 6dB ge-genuber dem Winkel θ = 111.6o und ist damit viermal so laut wie senkrecht zur Anstromnug.Der Pegel der Grundfrequenz sinkt um jeweils zwei dB zur vorigen Position. Die dritte Har-monische sinkt ebenfalls im Pegel. Bei einem Winkel θ = 158.4o ist der Absolutpegel vonGrundfrequenz und zweiter Harmonischer gleich hoch.


0.309 0.31 0.311 0.312 0.313 0.314 0.315 0.316t [s]

-2

-1

0

1

2

p’(t

) [P

a]

0.309 0.31 0.311 0.312 0.313 0.314 0.315 0.316t[s]

-2

-1

0

1

2

p’(t

) [P

a]


180o


0o


Fur den Winkel von θ = 180o, d.h in Richtung, aus der der Zylinder angestromt wird, istdas Signal in Abbildung 39 dargestellt. Des weiteren ist das Signal bei θ = 0o dargestellt,bei dem der Schalldruckpegel SPL am geringsten ist. Der Signal schwingt mit doppelterGrundfrequenz. Die Amplitude ist bei θ = 0o kleiner. Das Verhalten der Schallabstrahlungvor und hinter dem Zylinder ist ahnlich, die Pegel sind hinter dem Zylinder grundsatzlichgeringer. Die Spektren sind in Abbildung 40 dargestellt.

200 400 800 2k1k 3kf [Hz]

0

20

40

60

80

100

SP

L [d

B, r

el. 2

*10-5

Pa]

200 400 800 1k 2k 3k400 800 2k1k 3kf [Hz]

0

20

40

60

80

100

SP

L [d

B, r

el. 2

*10-5

Pa]

180o


0o


Der Pegel der Grundfrequenz ist beim Beobachter in Richtung θ = 180o knapp unter 70dBgesunken, der Pegel der zweiten Harmonischen auf fast 88dB angestiegen und somit bedeu-tend lauter als der Grundton. Der Unterschied betragt 18dB. Die dritte Harmonische ist aufca. 62dB gesunken. Die Verhaltnisse der Frequenzen zueinander sind in der gegenuberliegen-den Richtung ( θ = 0o) gleich.

Die Ergebnisse fur unterschiedliche akustische Lange lz unterscheiden sich nicht. Bei Steige-rung von lz steigt proportional der Pegel. Die Richtcharakteristik besitzt die gleiche Form,


nur mit einem Offset im Pegel. Die Pegel der Peakfrequenzen und hoheren Harmonischensteigen proportional zur Steigerung der effektiven akustischen Lange.Eine Aussage uber den absoluten Pegel ist mit der Simulation nicht moglich. Die Simulationermoglicht jedoch Aussagen zum Pegel der Frequenzen zueinander.

6.2.5 FW-H auf den durchlassigen Kontrollflachen fi ausgewertet

Die Geschwindigkeitskomponenten und Drucke auf den durchlassigen Kontrollflachen werdenzur Zeit entlang einer Linie aus den Stromungsdaten extrahiert. Die Schnittpunkte dieserLinie mit den Volumenelementen der Stromungssimulation sind die Koordinaten der Ele-mente der Kontrollflache. Die Werte auf den Elementen werden auf Basis der Großen derumgebenden finiten Volumen interpoliert.Im Rahmen der Auswertung kommt man jedoch zu dem Schluß, daß eine Interpolationder Daten zu sehr großen Fehlern fuhrt, da die Schwankungen des Druckes und der Ge-schwindigkeiten in der Zeit und in den benachbarten Volumenzellen starken Schwankungenunterworfen sind.Die Fehler sind so groß, daß keine sinnvollen Auswertungen moglich sind. Casalino u.a. [2]stellen diesen Effekt in ihrer neuesten Veroffentlichung dar. Die Fehler durch Interpolationder Werte von dem CFD-Gitter auf eine Kontrollflache liegen bei Casalino im Bereich von15dB bis 20dB.Eine komplette Neuprogrammierung der Extraktionsroutinen in ELAN2 und neue Simula-tionen der Umstromung des Zylinders fur die Berechnung der Daten auf den Kontrollflachensprengen den zeitlichen Rahmen der Studienarbeit. Auf eine Auswertung auf den durchlassi-gen Kontrollflachen muß deshalb im Rahmen dieser Arbeit verzichtet werden.


7 NACA 0015

Die Simulation eines symmetrischen NACA0015-Flugelprofils erfolgt als zweiter Testfall.Shen und Sørensen [15] veroffentlichen die CAA-Simulation (Computational Aeroacoustics)eines NACA0015-Flugelprofils bei einem Anstellwinkel von α = 20o und einer Reynoldszahlvon Re = 1.5 · 106 und einer Machzahl von M = 0.2. Die Daten von Shen werden zumVergleich herangezogen.

7.1 Stromungssimulation

7.1.1 Konfiguration

Die Stromungssimulation des NACA0015 erfolgt bei Re = 1.5 · 106 und einem Anstellwinkelvon α = 20o. Die vollturbulente Simulation basiert auf dem LLR-k−ω−Modell. Die Konvek-tionsterme werden mit Hilfe des TVD-MUSCLE-Schemas diskretisiert. Der Turbulenzgradder Grundstromung ist auf Tu = 1.0% eingestellt.

Die Rechnung ist instationar und geht uber 8000 Zeitschritte mit der dimensionslosen Zeit-schrittweite ∆t = 0.005. Die Gesamtsimulationszeit betragt TA = 40. Die Zeitdiskretisierungist von erster Ordnung genau.

Fur die akaustische Auswertung, werden die Ergebnisse der Simulation auf einen realenFlugel skaliert:Die Lange des Flugelprofils betragt l = 0.3353m bei einer Anstromgeschwindigkeit vonu∞ = 68m

s. Dies entspricht einer Machzahl von M = 0.2. Die kinematische Viskositat betragt

ν = 1.52 · 10−5 m2

sund die Dichte % = 1.2 kg

m3 .Damit ergibt sich die mit der Lange l gebildete Reynoldszahl Re = 1.5 · 106.


7.1.2 Rechengitter

Das Gitter besteht aus 24240 Volumenelementen in zwei Blocken. Die Volumenelemente sindzur Flugeloberflache stark verdichtet. Das Gitter im Block um den Flugel besitzt C-Struktur,an den sich der zweite Block anschließt.Das Gitter erlaubt das Rechnen mit LOW-REYNOLDS-Randbedingung, d.h. die Grenz-schicht wird komplett bis in den linearen Bereich aufgelost. Die normierte Hohe y+ = yuτ

ν

der wandnachsten Zellen liegt im Bereich von y+ ≤ 1.0.

Abbildung 41 zeigt das Rechengitter in unmittelbarer Umgebung des Flugelprofils, mit derAnstromung von links im Winkel von α = 20o.

Abbildung 41: Verwendetes Rechengitter um das NACA0015 Profil


7.1.3 Randbedingungen

Auf der Flugeloberflache wird Wandhaftung angenommen. Das Fluid bewegt sich an derWand mit der gleichen Geschwindigkeit wie die Korperoberflache.Am Einstromrand werden die Geschwindigkeiten und Turbulenzgroßen vorgegeben, am Aus-fluß wird eine konvektive Randbedingung verwendet.

7.1.4 Auswertung

Die Umstromung des NACA0015 bei Re = 1.5 ·106 ist bedingt durch den großen Anstellwin-kel durch Ablosen der Stromung vor dem großten Querschnitt des Profils gekennzeichnet.Infolgedessen entsteht ein Ruckstromgebiet auf dem Flugel, das in der Große schwankt. Ander Hinterkante des Profils losen in regelmaßigem Abstand Wirbel ab.In Abbildung 42 ist der instationare Verlauf des Auftriebsbeiwertes cl dargestellt. uber derdimensionslosen Zeit dargestellt.

100 110 120 130 140 150dimensionslose Zeit [t . uo/l]

0.6

0.8

1

1.2

1.4

1.6

cl

cl-Verlauf im ausgewerteten Zeitintervall T

A = 40

Ende der Startrechnung

Abbildung 42: Instationarer Verlauf des Auftriebswertes cl des NACA0015-Profils bei Re = 1.5 ·106, Machzahl M = 0.2, Anstellwinkel α = 20o

Das Signal ist periodisch mit einer Oberschwingung,die bei Maximadurchgang des Signalszu erkennen ist. Die Lange des Signals, das fur die Akustiksimulation genutzt wird, betragtin dimensionsloser Zeit TA = 40. Dimensionsbehaftet ergibt sich ein Zeitintervall von TA =0.1972s.Auf Basis des Auftriebsbeiwertes ist die Strouhalzahl der Wirbelablosung bestimmt. In Ab-bildung 43 ist die Fouriertransformierte des Auftriebsbeiwertes dimensionslos dargestellt.Die Strouhalzahl der Wirbelablosung betragt St = 0.58. Dies entspricht einer Frequenz vonf = 117.6Hz. Shen gibt eine Strouhalzahl von StShen = 0.862 an, jedoch bestatigen DES-Simulationen den geringeren Wert der aktuellen URANS-Simmulation.Die gemittelten Werte fur den Auftriebsbeiwert cl und Widerstandsbeiwert cd liegen beicl = 1.032 bzw. cd = 0.365. Der berechnete mittlere Auftriebsbeiwert des NACA0015 bei


0 1 2 3 4 5 6St

0

0.2

0.4

0.6

0.8

1

norm

iert

es S

pekt

rum

des

Auf

trie

bsbe

iwer

tes

c l

Abbildung 43: Spektrum des Auftribsbeiwertes des NACA0015-Profils bei Re = 1.5 · 106, Mach-zahl M = 0.2, Anstellwinkel α = 20o

einem Anstellwinkel von α = 20o stimmt mit experimentellen Werten von cl = 0.995 gutuberein.

Die globalen Werte der URANS-Rechnungen stimmen gut mit den experimentellen Datenuberein. Die vorhergesagten Schwankungen des Auftriebswertes cl und die Strouhalzahl Stsind mangels Vergleichsdaten in ihrer quantitativen Aussage nicht bestatigt.


7.2 Akustiksimulation

7.2.1 CAA-Vergleichsdaten

Shen und Sørensens [15] Rechnungen basieren auf einer instationaren, inkompressiblen CFD-Rechnung, die um einen Akustikansatz erweitert wird. Die Simulation des NACA0015-Profilsbei einer Reynoldszahl von Re = 1.5 · 106 und einem Anstellwinkel von α = 20o fuhren Shenu.a. [15] auf einem 121 × 161 O-Gitter mit y+ ≤ 4.0 mit einem Eingleichungsmodell vonBaldwin-Barth durch. Die Akustikrechnung erfolgt auf einem 101× 81 O-Gitter, mit einemVergleich der Ergebnisse auf einem 201× 81 O-Gitter fur die Akustiksimulation.In Abbildung 44 ist die von Shen errechnete Richtcharakteristik dargestellt.

Abbildung 44: Richtcharakteristik des NACA0015 im Abstand von r = 4.02m nach Shen [15],Re = 1.5 · 106, Machzahl M = 0.2, Anstellwinkel α = 20o

Die Pegel liegen im Bereich von 91dB bis 110dB. In Bewegungsrichtung ist das Flugelprofillauter.Shen stellt noch Zeitschriebe des Schalldruckes in r = 4.02m Entfernung dar. Die Ergebnissesind beschrankt auf die Winkel θ = 0o, 90o und 270o. Fur den Winkel θ = 90o gibt Shenebenfalls den Zeitverlauf mit feinerer Gitterauflosung an. In Abbildung 45 ist der Zeitverlaufdes normierten Schalldruckes p′/p0 fur zwei verschiedene Gitter dargestellt.Die dominierende Frequenz liegt nach Shen im Bereich der doppelten Strouhalzahl. Die Un-terschiede der Ergebnisse sind bei den beiden Gittern groß. Der Pegel wachst mit Erhohungder Auflosung des Gitters.Die Schalldruckverlauf nach Shen [15] senkrecht zur Anstromung des Profils, berechnet aufdem groben Gitter, ist in Abbildung 46 dargestellt.


Abbildung 45: Zeitverlauf des normierten Schalldruckes (p− p0)/p0 des NACA0015 im Abstandvon r = 4.02m nach Shen, Re = 1.5 · 106, Machzahl M = 0.2, Anstellwinkel α = 20o; a) θ = 0o auf101× 81-Gitter und b) θ = 0o auf 201× 81-Gitter

Abbildung 46: Zeitverlauf des normierten Schalldruckes (p− p0)/p0 des NACA0015 im Abstandvon r = 4.02m nach Shen [15] in Richtung θ = 90o und θ = 270o, Re = 1.5 ·106, Machzahl M = 0.2,Anstellwinkel α = 20o; 101× 81-Gitter

In Richtung θ = 90o dominiert nach Shen wiederum die doppelte Strouhalzahl, in Richtungθ = 270o ist jedoch die Grundfrequenz am starksten im Pegel vertreten.

7.2.2 Konfiguration

In Abbildung 47 sind die Kontrollflachen fi dargestellt. Auf diesen Flachen sind die insta-tionaren Druckdaten und Geschwindigkeitskomponenten aus der Stromungssimulation ex-trahiert. Die Oberflache f1 ist identisch mit der Oberflache des Flugelprofils, hier sind nurdie Druckdaten vorhanden. Die Kontrollflache f2 ist durchlassig, es werden die Druckdatenund Geschwindigkeitskomponenten zur Auswertung benutzt. Die Flache ist der Profilformangepaßt und in unmittelbarer Nahe zum Flugel positioniert.Die Parameter fur die Berechnungen mit C2Noise sind durch die Stromungssimulation vorge-geben: Die Kontrollflache und der Beobachter bewegen sich im Bezugssystem mit v = 68m

s


-0.2 -0.1 0 0.1 0.2 0.3 0.4x[m]

-0.25

-0.2

-0.15

-0.1

-0.05

0

0.05

0.1

0.15

0.2

y [m

]

f1: undurchlässige Flügeloberfläche

f2: durchlässige Kontrollfläche

Abbildung 47: Kontrollflachen fi, auf denen die Druck- und Geschwindigkeitsdaten aus derStromungssimulation extrahiert werden

nach links. Der Schalldruckverlauf wird an den Messpositionen von Shen in r = 4.02msenkrecht zur Anstromung berechnet. Zusatlich werden auf einem Kreis mit dem Radiusr = 4.02m um den Flugel herum an 50 gleichmaßig verteilten Punkten der Schalldruckver-lauf berechnet, um die Richtcharakteristik der Schallabstrahlung darstellen zu konnen.

Die effektive akustische Ausdehnung in z-Richtung lz der Rechnung in C2Noise betragtca. dreimal die Profillange. Die effektive akustische Gesamtlange ist lz = 1.02m, mit einerAuflosung ∆z = 0.02m und Anzahl der Elemente in z-Richtung nz = 25.

7.2.3 Auswertung der FW-H-Simulation

Alle folgenden Auswertungen sind fur Beobachter im Abstand von r = 4.02m. Die Gradanga-ben in den Abbildungen beziehen sich auf den Winkel zwischen der Hauptstromungsrichtungund der Richtung zum Beobachter. Die Angaben θ = 90o und θ = 270o sind senkrecht zurAnstromung, θ = 180o ist entgegen der Anstromung des NACA0015.

7.2.4 FW-H auf der festen Korperoberflache f1 ausgewertet

Die Richtcharakteristik der Schallabstrahlung des NACA0015-Profils auf Basis des FW-Hist in Abbildung 48 dargestellt.Der Richtcharachteristik ist symmetrisch, mit einem Schalldruckpegel im Bereich von 91.5dB ≤SPL ≤ 110dB. Der minimale Pegel tritt im Winkel von θ = 340o auf, d.h. direkt in derVerlangerung der Skelettlinie des Profils. Der Pegel im Bereich hinter dem Flugel ist durch-weg geringer als vor diesem. In Richtung des Anstellwinkels, d.h. bei θ = 160o, verringertsich der Schalldruckpegel auf SPLθ=160o = 96.5dB, ist damit jedoch ca. 5dB lauter als hinterdiesem.Der Schalldruckpegel ist unter dem Flugel im Schnitt ca. 2dB lauter als auf der Oberseite.


120 110 100 90 90 100 110 120

SPL [dB, rel. 2*10-5

Pa]

120

110

100

90

90

100

110

120NACA0015:

- auf fester Oberflächeausgewertet

- TA = 0.2s

Abbildung 48: Directivity-Plot des NACA0015, Re = 1.5 · 106 Anstellwinkel α = 20o

Die symmetrische Richtcharakteristik korrespondiert mit der alleinigen Auswertung derOberflachenintegrale auf der undurchlassigen Oberflache. Ist die Relativgeschwindigkeit zwi-schen Beobachter und Quelle null, liefert nur der “loading-noise”-Term einen Beitrag zurLosung. Dieser ist durch eine dipolartige Richtcharakteristik gekennzeichnet.Wie stark der Einfluß der Quellen im Raum ist und wie sich die Richtcharakteristik durchdiese andern, muß durch Auswertung der Kontrollflache f2 geklart werden.


In Abbildung 49 ist der Zeitverlauf des Schalldruckes p′ an den zu Shen (vergleiche Abbildung46) korrespondierenden Beoachterpositionen bei θ = 90o und θ = 270o dargestellt. Dasdargestellte Zeitintervall von T = 0.1s entspricht einem dimensionslosen Zeitintervall vonT = 20.3 und ist damit identisch zu den Abbildungen von Shen (Abbildungen 45 und 46).

0.6 0.62 0.64 0.66 0.68 0.7t [s]

-15

-10

-5

0

5

10

15

p’(t

) [P

a]

0.6 0.62 0.64 0.66 0.68 0.7t[s]

-15

-10

-5

0

5

10

15

p’(t

) [P

a]

90o


270o

Abbildung 49: Schalldruckverlauf p’ des NACA0015 bei Re = 1.5 · 106, Anstellwinkel α = 20o

Der Schalldruckverlauf an den gegenuber liegenden Positionen ist gegenphasig, der Signal-verlauf ist an den beiden Positionen ahnlich. Die Amplituden in Richtung θ = 270o sindgroßer. Die Signale sind durch die Grundfrequenz dominiert. Eine ausgepragte Spitze imSignal ist auf hohere Harmonische mit starken Pegeln zuruckzufuhren.

0.6 0.62 0.64 0.66 0.68 0.7t [s]

-15

-10

-5

0

5

10

15

p’(t

) [P

a]

0.6 0.62 0.64 0.66 0.68 0.7t[s]

-15

-10

-5

0

5

10

15

p’(t

) [P

a]

0o


340o


Der Schalldruckverlauf hinter dem Flugel bei θ = 0o und an der Position mit dem Minimal-pegel θ = 340o ist in Abbildung 50 dargestellt. Das Signal bei θ = 0o verhalt sich ahnlichdem bei θ = 90o, 270, jedoch mit deutlich geringeren Amplituden. Sowohl die Amplitude derGrundfrequenz, wie auch die durch die hoheren Harmonischen verursachte Spitze im Signalwerden kleiner.Bei θ = 340o ist kein Uberschwingen im Signal vorhanden, jedoch besitzt der Schalldruck-verlauf “Sagezahnform”, die hoheren Harmonischen sind mit deutlichem Pegel vertreten.


In den Abbildungen 51 bis 53 ist der Schalldruckverlauf von θ = 160o, der Richtung mit demgeringsten Pegel vor dem Flugel, bis in Richtung des maximalen Pegels bei θ = 250o in sechsaufeinanderfolgenden Richtungen dargestellt.

0.6 0.62 0.64 0.66 0.68 0.7t [s]

-15

-10

-5

0

5

10

15

p’(t

) [P

a]

0.6 0.62 0.64 0.66 0.68 0.7t[s]

-15

-10

-5

0

5

10

15

p’(t

) [P

a]

160o


166o


0.6 0.62 0.64 0.66 0.68 0.7t [s]

-15

-10

-5

0

5

10

15

p’(t

) [P

a]

0.6 0.62 0.64 0.66 0.68 0.7t[s]

-15

-10

-5

0

5

10

15

p’(t

) [P

a]

180o


202o


Die starke Zunahme der Amplituden der Grundfrequenz mit zunehmendem Winkel ist guterkennbar. Die Spitze im Signalverlauf, verursacht durch die hoheren Harmonischen, wachstanalog zur Amplitude der Grundfrequenz, ist jedoch erst bei θ = 223o und 250o stark aus-gepragt.


0.6 0.62 0.64 0.66 0.68 0.7t [s]

-15

-10

-5

0

5

10

15

p’(t

) [P

a]

0.6 0.62 0.64 0.66 0.68 0.7t[s]

-15

-10

-5

0

5

10

15

p’(t

) [P

a]

223o


250o


In den Abbildungen 54 bis 56 sind die Spektren des Schallruckes an den Positionen vonθ = 160o bis θ = 270o dargestellt.

50 100 200 400 800 1k800f [Hz]

40

60

80

100

120

SP

L [d

B, r

el. 2

*10-5

Pa]

50 100 200 400 800 1kf [Hz]

40

60

80

100

120

SP

L [d

B, r

el. 2

*10-5

Pa]

160o


180o

Abbildung 54: Spektrum des NACA0015 bei Re = 1.5 · 106, Anstellwinkel α = 20o

Bei θ = 160o, der Richtung mit dem geringsten Schalldruckpegel vor dem Flugel, liegt derPegel der Grundfrequenz f = 116Hz bei SPL116Hz = 99dB. Die zweite Harmonische liegt15dB die dritte und vierte Harmonische 19dB unter dem Pegel der Grundfrequenz.In Richtung θ = 180o steigt der Pegel der Grundfrequenz und der zweiten Harmonischengeringfugig um ca. 2dB. Die dritte und vierte Harmonische steigen stark an. Der Pegelbei f3.H. = 345Hz ist ca. 7dB leiser als der Grundton und damit lauter als die zweiteHarmonische. Die vierte Harmonische bei f4.H = 460Hz ist mit SPL4.H = 90dB ebenfallslauter als die zweite Harmonische mit SPL2.H = 87dB. Noch hohere Frequenzen sind mitdeutlich geringeren Pegeln vertreten.


50 100 200 400 800 1k800f [Hz]

40

60

80

100

120S

PL

[dB

, rel

. 2*1

0-5 P

a]

50 100 200 400 800 1kf [Hz]

40

60

80

100

120

SP

L [d

B, r

el. 2

*10-5

Pa]

202o


223o


Die gleichen Verhaltnisse sind auch bei ansteigenden Beobachterwinkeln θ = 202o undθ = 223o zu beobachten. Die Pegel aller Frequenzen steigen mit zunehmenden Winkelngleichmassig um jeweils rund 3dB.

50 100 200 400 800 1k800f [Hz]

40

60

80

100

120

SP

L [d

B, r

el. 2

*10-5

Pa]

50 100 200 400 800 1kf [Hz]

40

60

80

100

120

SP

L [d

B, r

el. 2

*10-5

Pa]

250o


270o


Bei θ = 250o, in Abbildung 56, erreicht der Pegel der Grundfrequenz mit 110dB sein Maxi-mum. Die hoheren Harmonischen steigen ebenfalls nur noch gering. Die Frequenzen großer500Hz sind im Verhaltnis starker gestiegen, die funfte Harmonische bei f5.H. = 575Hz er-reicht mit einem Pegel von SPL5.H = 96dB erstmals das Niveau der zweiten Harmonischen.Senkrecht zur Anstromung, bei θ = 270o, sind die Pegel ahnlich denen in Richtung θ = 250o.


In Abbildung 57 sind die Spektren hinter dem Flugel bei θ = 0o und θ = 340o dargestellt.

50 100 200 400 800 1k800f [Hz]

40

60

80

100

120

SP

L [d

B, r

el. 2

*10-5

Pa]

50 100 200 400 800 1kf [Hz]

40

60

80

100

120

SP

L [d

B, r

el. 2

*10-5

Pa]

0o


340o


In der Richtung mit dem geringsten Schalldruckpegel θ = 340o hat die Grundfrequenz f =115Hz einen Pegel von SPL = 94dB. Die hoheren Harmonischen fallen in der dB-bewertetenSkala linear mit der Frequenz ab.

Bei θ = 0o sind die Pegel aller Frequenzen mit Ausnahme der zweiten Harmonischen um 8dBgestiegen, der lineare Abfall zu hoheren Frequenzen bleibt erhalten. Der Pegel der zweitenHarmonischen ist gegenuber θ = 340o um 3dB gesunken.

Die Berechnung der Richtcharakteristik des NACA0015 auf Basis des FW-H deckt sichgrundsatzlich mit der Berechnung von Shen auf Basis einer CFD/CAA-Rechnung. Die domi-nierenden Frequenzen sind bei Shen jedoch die doppelte Strouhalzahl, in den vorliegendenRechnungen liegen sie identisch zur Strouhalzahl der Wirbelablosung. Ein Grund fur dieausgepragten Pegel bei der doppelten Strouhalzahl gibt Shen nicht an.

Die Spektren an den verschiedenen Beobachterpositionen geben einen guten Einblick in dasVerhalten der hoheren Frequenzen bei verschiedenen Beobachterwinkeln. Erst mit den Ein-zelspektren wird das haufig gegensatzliche Verhalten von Grundfrequenz zu den hoherenHarmonischen darstellbar.

In den Auswertungen sind wiederum nur tonale Anteile in den Spektren enthalten, da dieEingangsdaten des FW-H auf zweidimensionalen Stromungssimulationen beruhen.

7.2.5 FW-H auf der durchlassigen Kontrollflache f2 ausgewertet

Die Interpolation der Stromungsdaten der Stromungssimulation auf der Kontrollflache f2

fuhrt zu großen Fehlern, die eine sinnvolle Auswertung der Daten verhindern. Da eine Neu-berechnung im zeitlichen Rahmen dieser Arbeit nicht moglich ist, muß auf die Auswertungder durchlassigen Kontrollflache f2 verzichtet werden.

8 ZUSAMMENFASSUNG 85

8 Zusammenfassung

Im Rahmen dieser Arbeit ist auf Basis von zweidimensionalen, instationaren URANS-Stro-mungssimulationen mit dem Programm ELAN2 und dem Programm C2Noise der Schall-druck im Fernfeld berechnet worden.

Die Simulation der Schallausbreitung erfolgt in C2Noise durch die Auswertung der instati-onaren Druck- und Geschwindigkeitskomponenten auf den Kontrollflachen mit dem Algorith-mus von Ffowcs-Williams und Hawkings (FW-H). Der FW-H ist die Losung der Lighthill-Differentialgleichung im Fernfeld, die die Schallerzeugung und -ausbreitung von bewegtenKorpern in bewegten Medien beschreibt. Die Kontrollflache kann dabei mit der Oberflachedes umstromten Korpers ubereinstimmen oder diesen komplett einschließen. Alle akustischenQuellen innerhalb der Kontrollflache werden erfasst.

Die zweidimensionale Implementierung des FW-H basiert auf einem vorhandenem 3D-FW-H-Code. Bei Nutzung des dreidimensionalen FW-H, auf der Basis von Daten aus zweidimen-sionalen Stromungssimulationen, mußten im kompletten Programmcode die Schnittstellenund die Routinen zur Berechnung der Oberflachendaten modifiziert werden. Die Berechnungder Quellterme und die Bestimmung der Quellzeiten sind identisch zum 3D-Code. Bei derBerechnung der akustischen Quellen werden die Werte aus der zweidimensionalen Stromungs-simulation den Elementen im R3 zugeordnet. Die Extraktion der Stromungsdaten erfordertedie Modifikation von Routinen in ELAN2 und die Neuentwickelung von Tools zur Kon-vertierung und Dimensionalisierung der Stromungsdaten, in dem fur den FW-H optimalen,unformatierten Datenformat.

Die Effekte einer zweidimensionalen Quelle, d.h. in der z-Richtung identisch fortgesetztenQuelle, wurden an einem δ-Puls untersucht. Durch Variation der Ausdehnung der 2D-Quellein z-Richtung konnten die Haupteffekte von zweidimensionalen Quellen simuliert werden.Mit zunehmender Ausdehnung der Quelle steigt der Pegel der 2D-Quelle kontinuierlich an.Wird die Ausdehnung der Quelle weiter gesteigert, bildet sich eine Signalschleppe. Ein sehrkurzer, zweidimensionaler Puls erzeugt dann ein langes Zeitsignal.

Fur die Simulation des Schalls von realen umstromten Korpern ergeben sich folgende Ein-schrankungen: Die effektive akustische Lange lz, die fur die Simulation in C2Noise eingestelltwird, muß kleiner sein als die reale Ausdehnung des umstromten Korpers in z-Richtung. DiePegel werden sonst zu hoch vorhergesagt, da in z-Richtung die Phasenlage von Druck- undGeschwindigkeitskomponenten gleich ist. Bei real umstromten Korpern ist dies nicht derFall. Die zeitliche Auflosung der Simulation ist eingeschrankt durch die Diskretisierung derKontrollflache: ∆xi ≤ c · ∆t

2.

Der erste angewandte reale Testfall ist der umstromte Kreiszylinder. Die Stromung ist durchwechselseitiges gleichmaßiges Ablosen von Wirbeln gekennzeichnet. Die Stromungssimulati-on auf Basis von RANS-Rechnungen ist jedoch nur ungenugend in der Lage, die korrekteinstationare Ablosung im subkritischen Reynoldszahlbereich von stumpfen Korpern zu mo-dellieren. Die berechneten Frequenzen im Fernfeld auf Basis des FW-H entsprechen denender Stromungssimulation. Ist die instationare Stromungssimulation fehlerbehaftet, wird derFehler durch den FW-H auch am Beobachterstandort im Fernfeld zu beobachten sein. Ausdiesem Grund sind die simulierten Frequenzen im Vergeich zu den Experimenten zu hoch.

Die Ergebnisse auf Basis des quasi-zweidimensionalen FW-H sind stark tonal. Die breitban-

86 8 ZUSAMMENFASSUNG

digen Signale der Experimente, die auf dreidimensionalen Effekten beruhen, kann der 2D-Ansatz prinzipbedingt nicht wiedergeben. Der zweidimensionale Ansatz ermoglicht Aussagenuber den Einfluß der hoheren Harmonischen. Die Anderungen im Pegel der Oberwellen, jenach Beobachterposition, sind meist starker als die Anderungen im Pegel der Grundfrequenz.In experimentellen Untersuchungen sind die Auswirkungen der hoheren Harmonischen meistnicht eindeutig erkennbar, da diese haufig in den breitbandigen Signalen untergehen.Die Verhaltnisse der simulierten Pegel der Oberwellen zu der Grundfrequenz sind mit demverfolgten Ansatz in guter Ubereinstimmung mit den Experimenten. Hier liegt die Haupt-anwendung des zweidimensionalen Ansatzes.Die Simulation des Schalls an beliebigen Beobachterpukten, die in Experimenten nur schwie-rig realisierbar sind, ermoglicht beispielsweise die Berechnung der Richtcharakteristik derSchallabstrahlung des umstromten Zylinders. Dabei ist es auch moglich, die Pegel der ver-schiedenen Frequenzen in Abhangigkeit vom Beobachterstandort zu bestimmen. Das Verhal-ten der berechneten Pegel von Grundfrequenz zur ersten Oberschwingung ist im Falle desZylinders gegenlaufig.Der zweite Testfall ist die Simulation eines NACA0015 Profils. Dieser bestatigte die beim Zy-linder gewonnenen Erkenntnisse. Die Simulationen des Flugelprofils zeigt den starken Einflußder hoheren harmonischen Frequenzen auf die Losung. Die Richtcharakteristiken auf Basisder instationaren Druckdaten auf der Korperoberflache sind bezogen auf den Gesamtschall-druckpegel nahezu symmetrisch und sind dabei zusatzlich von der Frequenz abhangig.Der angestrebte Vergleich von Simulationen auf den festen Korperoberflachen mit frei de-finierten Kontrollflachen zur Erfassung der raumlichen Quellen konnte nicht durchgefuhrtwerden. Durch die Interpolation der Stromungsdaten auf die Kontrollflachen entstehen sehrgroße Fehler, die eine Auswertung verhindern.

Die Simulation der Schallabstrahlung auf Basis von zweidimensionalen Stromungsdaten istmoglich. Es sind Aussagen uber das tonale Verhalten und das Verhaltnis der Frequenzen zu-einander moglich. Die Qualitat der Stromungssimulation hat entscheidenden Einfluß auf densimulierten Schall im Fernfeld. Die berechneten Frequenzen korrelieren bei den berechneteneinfachen Stromungskonfigurationen mit den Schwankungen des Auftriebsbeiwertes.Die Verringerung des numerischen Aufwandes durch den zweidimensionalen Ansatz sindgroß. Die instationaren Stromungssimulationen von ca. 10000 Zeitschritten sind bei rea-litatsnahen Reynoldszahlen auf einem PC durchgefuhrt worden. Die Rechenzeiten lagen imBereich von einer bis maximal zwei Wochen bei den Teilrechnungen. Die Berechnung desSchalls im Fernfeld erfordert pro Beobachterposition zwischen einer und funf Minuten.Auf Basis dieser Erkenntnisse sind weitere Arbeiten zur Validierung des Einsatzes des zwei-dimensionalen FW-H bei realen Stromungskonfigurationen notig. Die Berechnung auf durch-lassigen Kontrollflachen erfordert eine komplett neue Implementierung der Extraktionsrou-tionen in ELAN2 auf direkter Basis der finiten Volumen ohne Interpolation. Nur auf diesemWege ist die Untersuchung von komplizierteren Stromungskonfigurationen moglich.Ein Vergleich mit dreidimensionalen CFD/FW-H -Rechnungen und experimentellen Datenist zur Bestimmung der Aussagekraft von zweidimensionalen Akustiksimulationen sinnvoll.Der Einfluß der raumlichen Quellen und damit der Einfluß der Lage der Kontrollflachen aufdie Ergebnisse sollte untersucht werden.

8 ZUSAMMENFASSUNG 87

Literatur

[1] K.S. Brentner and F. Farassat. Analytical comparison of the acoustic analogy andkirchhoff formulation for moving surfaces. AIAA Journal, 36(8):1379–1386, 1998.

[2] D Casalino, M. Jacob, and M. Roger. Prediction of rod-airfoil interaction noise usingthe FW-H analogy. AIAA Journal, (2543), 2002.

[3] J.S. Cox, C.L. Brentner, K.S.and Rumsey, and B.A. Younis. Computation of vor-tex shedding and radiated sound for a circular cylinder. 4th International Symposiumon Fluid-Structure Interactions, Aeroelasticity, Flow-Induced Vibration and Noise, AD-53(1):447–454, 1997.

[4] K. Ehrenfried. Skript - Stromungsakustik. http://vento.pi.tu-berlin.de, 20.5.2002.

[5] J.E. Ffowcs Williams and D.L. Hawkings. Sound generated by turbulence and surfacesin arbitrary motion. Philosophical Transactions of the Royal Society, A264:321–342,1969.

[6] W.F. King III, E. Pfizenmaier, M. Herrmann, and Neuhaus L. An experimental stu-dy of sound generated by flow interactions with cylinders. internal report, DeutscheForschungsanstalt fur Luft- und Raumfahrt e.V. (DLR), 1996.

[7] M.J. Lighthill. On sound generated aerodynamically, i: General theory. Proceedings ofthe Royal Society London, A221:564–587, 1952.

[8] D.P. Lockard. An efficient, two-dimensional implementation of Ffowcs Williams andHawkings Equation. Journal of Sound and Vibration, 229:897–911, 2000.

[9] T. Rung. Realizability linearer Stress-Strain Beziehungen. Institutsbericht, Nr. 04/98,Hermann-Fottinger-Institut, 1998.

[10] T. Rung. Entwicklung anisotroper Wirbelzahigkeitsbeziehungen mit Hilfe von Projekti-onstechniken. PhD thesis, Technische Universitat Berlin, 2000.

[11] T. Rung, D. Eschricht, J. Yan, and F. Thiele. Sound radiation of the vortex flow pasta generic side mirror. American Institute of Aeronautics and Astronautics Paper, 2340,2002.

[12] T. Rung and F. Thiele. Computational modelling of complex boundary-layer flows.In 9th Int. Symp. on Transport Phenomena in Thermal-Fluid Engineering, Singapore,1996.

[13] G. Schewe. On the force fluctuations acting on a circular cylinder in cross-flow fromsubcritical to transcritical reynolds numbers. Journal of Fluid Mechanics, 133:265–285,1983.

[14] S. Schmidt. Grobstruktursimulation turbulenter Stromungen in komplexen Geometrienund bei hohen Reynoldszahlen. Mensch und Buch Verlag, Berlin, 2002.

88 8 ZUSAMMENFASSUNG

[15] W.Z. Shen and J.N. Sorensen. Aeroacoustic modeling of turbulent airfoil flows. AIAAJournal, 39(6):1057–1064, 2001.

Berechnung der Schallabstrahlung umstr¨omter K ¨orper auf ...

Documents

Transcript of Berechnung der Schallabstrahlung umstr¨omter K ¨orper auf ...