Lineare Stabilitätstheorie

Die Lineare Stabilitätstheorie (kurz LST) beschreibt physikalisch in einer Strömung das Anwachsen wellenförmiger Störungen mit kleiner Amplitude. Durch Vorgabe eines stationären Strömungsfeldes lassen sich Anfachungsraten und Form der linearen Störungen in Abhängigkeit von Frequenz oder Wellenzahl bestimmen. Damit lässt sich gut erklären, warum sich für manche Strömungen kein stationärer Zustand einstellt. Trotz des lokalen Ansatzes kann mit ihr der anfängliche Bereich des laminar-turbulenten Umschlags bereits gut beschrieben werden. Die e^N-Methode, die eine gängige Methode zur Abschätzung der Position des laminar-turbulenten Umschlags ist, basiert auf den Anfachungsraten aus der linearen Stabilitätstheorie.

Geschichte

Viele in Natur und Technik auftretende Strömungen sind turbulent. Nur bei kleiner Reynolds-Zahl ergibt sich eine laminare Strömung. Lange Zeit war unklar, warum eine laminare Strömung umschlägt und sich stromab ein turbulenter Zustand einstellt. Eine der ersten Hypothesen war bereits, dass die laminare Strömung instabil sei. Ab den 1880er Jahren veröffentlichte John William Strutt, 3. Baron Rayleigh, verschiedene Ergebnisse zur Instabilität reibungsfreier Strömungen und entdeckte, dass ein Strömungsfeld beim Auftreten eines Wendepunkts im Geschwindigkeitsprofil instabil wird. Damit konnte zwar bereits das instationäre Verhalten eines Freistrahls erklärt werden, jedoch blieb die Fragestellung, warum eine Grenzschicht turbulent wird, weiter offen. Geoffrey Ingram Taylor hatte bereits 1915 den Verdacht, dass die Reibung auch destabilisierend wirken könnte. Typischerweise wurde jedoch von einem stabilisierenden Einfluss der viskosen Terme ausgegangen, sodass Fritz Noether 1921 aufgrund der bisherigen Forschungsergebnisse den Ansatz kleiner Störungen als nicht erfolgreich bewertete. Bemerkenswerterweise stellte im gleichen Jahr Ludwig Prandtl ähnliche Vermutungen wie Taylor an und initiierte die Forschung zur reibungsbehafteten Stabilitätstheorie in Göttingen. Werner Heisenberg beschäftigte sich im Rahmen seiner Dissertation mit dem laminar-turbulenten Umschlag, aber erst Walter Tollmien gelang es, die reibungsbehaftete Stabilitätstheorie zu entwickeln. Eine der ersten Anwendungen der linearen Stabilitätstheorie erfolgte 1933 durch Hermann Schlichting. Vor Beginn des Computerzeitalters führte Pretsch 1942 die ersten systematischen Stabilitätsrechnungen für Grenzschichten der Falkner-Skan Familie durch.

Außerhalb Deutschlands fand die lineare Stabilitätstheorie anfangs wenig Beachtung, da die vorausgesagten Instabilitätswellen nicht gemessen werden konnten. Außerdem überwog Skepsis darüber, dass ein linearer Ansatz sinnvolle Aussagen für ein hochgradig nichtlineares Phänomen wie die Turbulenz ermöglichen sollte. Erst durch den experimentellen Nachweis von Instabilitätswellen in einer Plattengrenzschicht durch Schubauer und Skramstadt Anfang der 1940er wurde der Göttinger Ansatz bestätigt (aufgrund der kriegsbedingten Zensur erfolgte eine Veröffentlichung der Ergebnisse 1947). Smith und Gamberoni sowie Van Ingen entwickelten aus der linearen Stabilitätstheorie die e^N-Methode, die auch heutzutage noch zur Transitionsvorhersage eingesetzt wird. Mit Beginn des Computerzeitalters in den 1960er Jahren wurden die Stabilitätseigenschaften zahlreicher Grenzschichtkonfigurationen berechnet. In Überschallgrenzschichten konnte Mack 1982 zusätzliche Instabilitätsmoden mittels der linearen Stabilitätstheorie nachweisen, die ihm zu Ehren auch als Mack-Moden bezeichnet werden.

Grundlagen der linearen Stabilitätstheorie

Die lineare Stabilitätstheorie betrachtet die Stabilität einer Strömung gegenüber kleinen Störungen. Hierbei wird davon ausgegangen, dass das Strömungsfeld in Querrichtung konstant ist. Für die Grenzschicht bedeutet dies, dass an einer ebenen Platte mit der Strömungsrichtung und der wandnormalen Richtung die Grundströmung konstant über der spannweitigen Richtung ist, was einer unendlichen Ausdehnung in -Richtung entspricht. Um das Verhalten kleiner Störungen mathematisch erfassen zu können, werden die Strömungsgrößen in den Navier-Stokes-Gleichungen in eine vorzugebende stationäre Grundströmung und instationäre Störgrößen unterteilt. Für die drei Geschwindigkeitskomponenten $u,v,w$ in -, - und -Richtung sowie Dichte $\rho$ , Druck und Temperatur ergibt sich somit:

$u=U_{0}(x,y)+u'(x,y,z,t),\quad v=V_{0}(x,y)+v'(x,y,z,t),\quad w=W_{0}(x,y)+w'(x,y,z,t)$

$\rho =\rho _{0}(x,y)+\rho '(x,y,z,t),\quad p=P_{0}(x,y)+p'(x,y,z,t),\quad T=T_{0}(x,y)+T'(x,y,z,t)$

Annahmen

Darauf aufbauend werden folgende Annahmen getroffen:

Die stationäre Grundströmung erfüllt die Navier-Stokes-Gleichungen. Dadurch verschwinden in den Navier-Stokes-Gleichungen alle Terme, die keine Störgrößen enthalten.
Parallele Grundströmung. Dabei wird z.B. das Aufdicken einer Grenzschicht vernachlässigt. Aus der Kontinuitätsgleichung folgt daraus, dass für die Grundströmung die normale Geschwindigkeitskomponente V₀ vernachlässigt wird. Dadurch ergibt sich eine lokale Theorie, bei der jede Position x in Stromabrichtung separat betrachtet wird. Daraus folgt: $U_{0}=U_{0}(y),\quad W_{0}=W_{0}(y),\quad P_{0}=P_{0}(y),\dots \quad V_{0}=0$
Kleine Störungen. Unter der Annahme kleiner Störungen können die Gleichungen bezüglich der Grundströmung linearisiert werden. Damit verschwinden alle Quadrate der Störgrößen aus den Gleichungen.

Störansatz

Für die Störgrößen wird ein Wellenansatz angesetzt, der z.B. für , die Störgeschwindigkeit in -Richtung, folgende Form hat:

$u'={\hat {u}}(y)\cdot e^{{i(\alpha x+\gamma z-\omega t)}}$

Dies entspricht einer Welle mit den Wellenzahlen $\alpha _{r}$ und $\gamma$ in - beziehungsweise -Richtung mit der Frequenz $\omega _{r}$ , deren Amplituden- und Phasenverlauf eine Funktion von sind. Die einzelnen Größen sind im Allgemeinen komplex und haben folgende physikalische Bedeutung:

$\alpha _{r}$ = Wellenzahl in Strömungsrichtung (Wellenlänge $\lambda _{x}=2\pi /\alpha _{r}$ )
$\gamma$ = Wellenzahl in Querrichtung (Wellenlänge $\lambda _{z}=2\pi /\gamma$ )
$k={\sqrt {\alpha _{r}^{2}+\gamma _{r}^{2}}}$ = Wellenzahl in Ausbreitungsrichtung
$\omega _{r}$ = Kreisfrequenz (Periodendauer $\lambda _{T}=2\pi /\omega _{r}$ )
$\alpha _{i}$ = räumliche Anfachungsrate in -Richtung (Anfachung für $\alpha _{i}$ < 0)
$\omega_i$ = zeitliche Anfachungsrate (Anfachung für $\omega_i$ > 0)
$\left|{\hat {u}}(y)\right|$ = Amplitudenverlauf von über
$\arg\{{\hat {u}}(y)\}$ = Phasenverlauf von über
$\omega _{r}/\alpha _{r}$ = Phasengeschwindigkeit der Störwelle in -Richtung
$\partial \omega _{r}/\partial \alpha _{r}$ = Gruppengeschwindigkeit der Störwelle in -Richtung

Die spannweitige Wellenzahl $\gamma$ ist aufgrund der unendlichen Ausdehnung der Strömung in -Richtung eine reelle Größe, da ein Imaginärteil $\gamma_i$ einem räumlichen Wachstum in spannweitiger Richtung entspräche.

Zeitliches und räumliches Modell

Bei der linearen Stabilitätstheorie wird zwischen zeitlichem und räumlichem Modell unterschieden. Beim zeitlichen Modell werden die reelle Wellenzahlen $\alpha _{r}$ und $\gamma$ vorgegeben, woraus sich die Frequenz $\omega _{r}$ und die zeitliche Anfachungsrate $\omega_i$ ergibt. Beim räumlichen Problem erhält man durch Vorgabe der Querwellenzahl $\gamma$ und der reellen Frequenz $\omega _{r}$ eine die Wellenzahl in -Richtung $\alpha _{r}$ und die räumliche Anfachungsrate $\alpha _{i}$ .

Ob das zeitliche oder das räumliche Problem anzuwenden ist hängt von der jeweiligen Problemstellung ab: handelt es sich um eine konvektive Instabilität, ist das räumliche Problem zu wählen, bei einer absoluten Instabilität das zeitliche. In vielen Fällen ist die räumliche Betrachtungsweise realistischer, da z.B. in einer Grenzschicht die Störungen während des Anwachsens stromab konvektiert werden, was einem Anwachsen im Raum entspricht.

Gaster-Transformation

Die Gaster-Transformation, ermöglicht die Transformation zwischen zeitlichem und räumlichem Modell für nahezu neutrale Wellen, das heißt unter der Annahme kleiner Anfachungsraten:

${\frac {\omega _{i}^{{\mathrm {zeit}}}}{\alpha _{i}^{{\mathrm {raum}}}}}=-{\frac {\partial \omega _{r}}{\partial \alpha _{r}}}$

wobei angenommen wird, dass sich Frequenz und -Wellenzahl im räumlichen und zeitlichen Modell nicht wesentlich unterscheiden:

$\alpha _{r}^{\mathrm {raum} }\approx \alpha _{r}^{\mathrm {zeit} },\quad \omega _{r}^{\mathrm {raum} }=\omega _{r}^{\mathrm {zeit} }$

Inkompressible Stabilitätsgleichungen

Aufgrund der elliptischen Form der Kontinuitätsgleichung für inkompressible Strömungen, ergeben sich bei Vernachlässigung von Kompressibilitätseffekten die Orr-Sommerfeld- und Squire-Gleichung zur Beschreibung der linearen Stabilität. Diese bilden zusammen ein Gleichungssystem sechster Ordnung.

Orr-Sommerfeld-Gleichung

Für inkompressible Strömungen wird das lineare Stabilitätsproblem durch die Orr-Sommerfeld-Gleichung beschrieben (nach William McFadden Orr 1907, Arnold Sommerfeld 1908):

${\begin{aligned}0=&\left(\alpha u_{0}+\gamma w_{0}-\omega \right)\left({\hat {v}}_{{yy}}-\left(\alpha ^{2}+\gamma ^{2}\right){\hat {v}}\right)\,-\,\left(\alpha u_{{0yy}}+\gamma w_{{0yy}}\right){\hat {v}}\\&{}+{\frac {i}{Re}}\left({\hat {v}}_{{yyyy}}-2\left(\alpha ^{2}+\gamma ^{2}\right){\hat {v}}_{{yy}}+\left(\alpha ^{2}+\gamma ^{2}\right)^{2}{\hat {v}}\right)\end{aligned}}$

Bei der Orr-Sommerfeldgleichung handelt es sich um eine gewöhnliche Differentialgleichung vierter Ordnung mit der komplexen Amplitudenfunktion ${\hat {v}}$ als einzige Variable. Man kann die Gleichung auch als Eigenwertproblem mit dem Eigenwert $\omega$ und dem Eigenvektor ${\hat {v}}(y)$ verstehen. Da bis zu vierte Ableitungen der zu lösenden Größe auftreten, sind vier Randbedingungen erforderlich.

Handelt es sich bei dem Rand um eine Wand, so ergeben sich aus der Haftbedingung und der Kontinuitätsgleichung:

${\hat {v}}(y_{{\mathrm {Wand}}})=0,\quad {\hat {v}}_{y}(y_{{\mathrm {Wand}}})=0$

Im Fall eines Freistromrandes muss die komplexe Amplitudenfunktion zum Fernfeld hin abklingen. Man kann zeigen, dass die Abklingrate in -Richtung der Wellenzahl in Stromabrichtung entspricht.

${\hat {v}}(y=\infty )=0,\quad {\hat {v}}_{y}(y=\infty )=0$

Squire-Gleichung

Die Lösung der Orr-Sommerfeld-Gleichung liefert für die Störung der normalen Geschwindigkeitskomponente Amplituden- und Phasenverlauf. Für einen rein zweidimensionalen Fall ( W_0=0 , $\gamma =0$ ) lassen sich die Verläufe der anderen Geschwindigkeitskomponente aus der Kontinuitätsgleichung berechnen. Im allgemeinen dreidimensionalen Fall ist jedoch die Squire-Gleichung zur Berechnung der Querströmungskomponenten erforderlich:

$i\left(\alpha u_{0}+\gamma w_{0}-\omega \right){\hat {\Omega }}-{\frac {1}{Re}}\left({\hat {\Omega }}_{{yy}}-\left(\alpha ^{2}+\gamma ^{2}\right){\hat {\Omega }}\right)=i\left(\alpha w_{{0y}}-\gamma u_{{0y}}\right){\hat {v}}$

mit der (wand)normalen Wirbelstärke

${\hat {\Omega }}=i(\gamma {\hat {u}}-\alpha {\hat {w}})$

als zu lösender Größe aus der sich die Störgrößen der Geschwindigkeitskomponenten und berechnen lassen. Hierbei ist der sich aus der Orr-Sommerfeld-Gleichung ergebende Verlauf von ${\hat {v}}$ auf der Rechthandseite zu verwenden. Die homogenen Lösungen der Squire-Gleichung (mit ${\hat {v}}=0$ ) sind ebenfalls gültige und physikalisch sinnvolle Lösungen und werden Squire-Moden genannt. Squire wies nach, dass diese Moden stets gedämpft sind.

Kompressible Stabilitätsgleichungen

Berücksichtigt man die Kompressibilität, so lässt sich mit der Kontinuitätsgleichung nicht eine Größe aus den Gleichungen eliminieren, wie dies bei der inkompressiblen Betrachtungsweise der Fall war. Mit den oben genannten Annahmen sowie dem Störansatz ergibt sich ein Gleichungssystem achter Ordnung, welches aus fünf Gleichungen besteht. Entsprechend der Ordnung des Gleichungssystems sind insgesamt acht Randbedingungen zu erfüllen:

Für die Wandrandbedingung gilt wie im inkompressiblen Fall die Haftbedingung. Unter der Annahme, dass die Wand thermisch träge ist und Temperaturfluktuationen mit Frequenzen von bis zu vielen Kilohertz, wie sie z.B. in Überschallgrenzschichten auftreten, nicht folgen kann, wird die Störgröße der Temperatur an der Wand ebenfalls zu null gesetzt.

${\hat {u}}(y_{{\mathrm {Wand}}})=0,\quad {\hat {v}}(y_{{\mathrm {Wand}}})=0,\quad {\hat {w}}(y_{{\mathrm {Wand}}})=0,\quad {\hat {T}}(y_{{\mathrm {Wand}}})=0$

Am Freistromrand gilt ähnlich wie im inkompressiblen Fall, dass die Störungen nach außen hin abklingen müssen. Allerdings lässt sich eine direkte Korrelation zwischen Abklingrate und Wellenzahl nicht aufstellen. Ist der Freistromrand weit genug entfernt, so lassen sich in guter Näherung die Störgrößen auch direkt zu null setzen.

Erkenntnisse aus der linearen Stabilitätstheorie

Aus der linearen Stabilitätstheorie lassen sich wesentliche physikalische Erkenntnisse über die Stabilität von Strömungen und damit über die Ursache des laminar-turbulenten Umschlags ableiten.

Reibungsfreie Instabilität

Für reibungsfreie Strömungen ( $Re\rightarrow \infty$ ) lassen sich aus der linearen Stabilitätstheorie folgende Theoreme ableiten:

Rayleigh-Theorem Nr. 1:

Eine notwendige Bedingung für reibungsfreie Instabilität ist ein Wendepunkt im Geschwindigkeitsprofil.

Rayleigh-Theorem Nr. 2:

Die Phasengeschwindigkeit einer angefachten Störung liegt stets zwischen dem Minimal- und dem Maximalwert der Grundströmung u(y).

Tollmien-Theorem:

Für eine Grenzschicht ist es notwendig und hinreichend für reibungsfreie Instabilität, dass die Grundströmung einen Wendepunkt besitzt.

Ein Wendepunkt beeinflusst somit wesentlich das Stabilitätsverhalten. Aus dem Theorem von Tollmien folgt außerdem, dass die Blasius-Grenzschicht erst durch den Einfluss der Reibung instabil wird.

Grenzschichtinstabilitäten

Stabilitätsdiagramm für die Blasius-Grenzschicht.

Führt man eine Stabilitätsanalyse für verschiedene Frequenzen an mehreren Positionen in Stromabrichtung durch, so erhält man ein Stabilitätsdiagramm. Für die inkompressible Blasius-Grenzschicht ist ein solches Stabilitätsdiagramm für zweidimensionale Wellen ( $\gamma =0$ ) rechts dargestellt. Dabei ist der Bereich der Anfachung in Abhängigkeit von der Frequenz $\omega$ und von der Reynolds-Zahl, gebildet mit der lokalen Verdrängungsdicke $\delta_1$ , gegeben. Die Neutralkurve (hier blau) gibt den Bereich an, in dem Störungen gedämpft beziehungsweise angefacht sind. Man sieht, dass ab einer Reynolds-Zahl $Re_{\delta _{1}}=520$ , der kritischen Reynolds-Zahl, erstmals Anfachung von Störungen existiert. Weiter stromab verschiebt sich der Bereich der Anfachung zu niedrigeren Frequenzen. Aufgrund der Form der Neutralkurve im Stabilitätsdiagramm ist oft auch von der Stabilitätsbanane die Rede.

Spektrum der Blasius Grenzschicht.

Das zeitliche Spektrum für die Blasiusgrenzschicht an der Stelle $Re_{\delta _{1}}=1000$ ist in der rechten Abbildung für eine Wellenzahl $\alpha _{r}=0{,}3$ dargestellt. Es existiert ein einzelner angefachter Eigenwert ( $\omega_i$ ) bei einer Frequenz $\omega _{r}=0{,}135$ , der der Tollmien-Schlichting-Welle (kurz TS-Welle) entspricht. Rechts des angefachten Eigenwertes befindet sich das kontinuierliche Spektrum, dessen Störungen alle gedämpft sind. Das kontinuierliche Spektrum entspricht Störungen in der Potentialströmung, entsprechend ist ihre Phasengeschwindigkeit gleich der dimensionslosen Freistromgeschwindigkeit ( $U_{\infty }=1$ ). Unterhalb des angefachten Eigenwertes befinden sich die gedämpften Squire-Moden, die sich aus der homogenen Lösung der Squire-Gleichung ( ${\hat {v}}=0$ ) ergeben und somit nur Störungen in und aufweisen.

Eigenfunktion des angefachten Eigenwertes: U-Profil der Tollmien-Schlichting-Welle.

Die zum angefachten Eigenwert gehörende Eigenfunktion zeigt für die Amplitude der Geschwindigkeitskomponente den typischen Verlauf einer TS-Welle mit dem globalen Maximum innerhalb der Grenzschicht und einem weiteren lokalen Maximum am Grenzschichtrand. Nach außen hin klingt die Störung exponentiell ab. Im Bereich des Nulldurchgangs der Amplitude existiert ein Phasensprung um $\pi$ .

Einfluss des Schräglaufwinkels

Der Schräglaufwinkel $\psi =\operatorname {atan}(\gamma /\alpha _{r})$ gibt den Winkel zwischen Ausbreitungsrichtung der Störwelle und der Langsströmungsrichtung an. Für den Einfluss des Schräglaufwinkels auf die Anfachungsraten in einer inkompressiblen Strömung gilt:

Squire-Theorem:

Für eine inkompressible zweidimensionale Grundströmung ergibt sich die kleinste kritische Reynoldszahl für zweidimensionale Störwellen.

Einfluss des Druckgradienten

Aus der linearen Stabilitätstheorie folgt, dass ein negativer Druckgradient $\partial p/\partial x<0$ die Strömung stabilisiert, ein Druckanstieg in Strömungsrichtung dagegen destabilisierend wirkt. Die Ursache hierfür liegt an der Form der Geschwindigkeitsprofile: Ein negativer Druckgradient bewirkt ein fülligeres Geschwindigkeitsprofil, wohingegen ein starker positiver Druckgradient einen Wendepunkt im u-Profil verursacht.

Überschallgrenzschichten

Skizze des lokalen Überschallgebiets der Phasengeschwindigkeit.

Mit zunehmender Strömungsgeschwindigkeit spielen Kompressibilitätseffekte zunehmend eine Rolle. Ab einer Machzahl von Ma = 3 existiert für die über der wandnormalen Richtung konstante Phasengeschwindigkeit an der Wand ein lokales Überschallgebiet, in dem die Phasengeschwindigkeit größer als die Summe aus lokaler Strömungs- und Schallgeschwindigkeit ist. Dies bewirkt das Auftreten einer zusätzlichen Instabilität, die üblicherweise als zweite Mode bezeichnet wird. Kennzeichnend für die zweite Mode ist eine Nullstelle in der Eigenfunktion des Druckes. Die Dicke des lokalen Überschallgebietes der Phasengeschwindigkeit ist hierbei ein Maß für die Anfachung der zweiten Mode.

Eine Erhöhung der Machzahl bewirkt in der Regel eine kleinere kritischere Reynoldszahl, deren Minimalwert bei einer Machzahl von ungefähr 4,5 erreicht wird. Ab einer Machzahl von 5 verbinden sich die Instabilitätsbereiche von erster und zweiter Mode im Stabilitätsdiagramm.

Die Kompressibilität hat auch Einfluss auf den Ausbreitungswinkel der maximal angefachten Störung. Während im inkompressiblen Fall für eine zweidimensionale Grenzschicht die angefachteste Störung einen Schräglaufwinkel von null Grad aufweist, nimmt er im Überschall von ca. 40° bei Mach 1 bis zu seinem Maximalwert von 65° bei Mach 3 zu.

Durch die Wandtemperatur kann das Stabilitätsverhalten der Grenzschicht stark beeinflusst werden: durch Kühlung ist es möglich, die erste Mode zu stabilisieren, allerdings führt dies zu größeren Anfachungsraten bei der zweiten Mode.

Numerische Lösungsverfahren

Für die Lösung der linearen Stabilitätsgleichungen gibt es zwei numerische Methoden, das Matrixverfahren und das Schießverfahren. Da die Grundströmungsprofile vorgegeben werden, können bei beiden Methoden die auftretenden y-Ableitungen der Grundströmung als bekannt angesehen werden, da sie z.B. mit Finiten Differenzen berechnet werden können.

Matrix-Verfahren

Beim Matrixverfahren, oft auch als direktes Verfahren bezeichnet, werden die zu lösenden Stabilitätsgleichungen in Matrixform dargestellt. Dabei werden die einzelnen Terme nach den -Ableitungen der zu lösenden Größen sortiert. Daraus ergibt sich ein Eigenwertproblem mit dem Eigenwert $\omega$ . Die Eigenvektoren bestehen aus den Werten der komplexen Amplitudenfunktionen an den diskreten -Stellen. Da die viskosen Terme der Navier-Stokes-Gleichungen auch zweite -Ableitungen enthalten, treten in der linearen Stabilitätstheorie Terme mit $\alpha ^{2}$ auf. Somit existiert das Eigenwertproblem nur für das zeitliche Modell, bei der räumlichen Betrachtungsweise wäre ein quadratisches Eigenwertproblem zu lösen.

Bei der zweidimensionalen Orr-Sommerfeld-Gleichung werden die Terme entsprechend den -Ableitungen von ${\hat {v}}$ und dem Auftreten des Eigenwertes $\omega$ sortiert:

$\left[-\alpha ^{3}U_{0}-\alpha U_{{0yy}}+{\frac {i\alpha ^{4}}{Re}}\right]\cdot \overline {v}+\left[\alpha U_{0}-{\frac {i\cdot 2\alpha ^{2}}{Re}}\right]\cdot \overline {v}_{{yy}}+\left[{\frac {i}{Re}}\right]\cdot \overline {v}_{{yyyy}}=\omega \left[1\right]\cdot \overline {v}_{{yy}}+\omega \left[-\alpha ^{2}\right]\cdot \overline {v}$

Hierbei steht $\overline{v}$ für den komplexen Eigenvektor, der aus den Werten von

${\hat {v}}(y)$ an den diskreten -Stellen besteht.

Die -Ableitungen von $\overline{v}$ lassen sich durch Finite Differenzen bestimmen. Somit ist die -Ableitung von $\overline{v}$ eine Linearkombination der einzelnen diskreten Werte, weshalb sich die -Ableitungen durch eine Matrixoperation ausdrücken lassen:

$\overline {v}_{{yy}}=\overline {D}_{2}\cdot \overline {v}$

$\overline {v}_{{yyyy}}=\overline {D}_{4}\cdot \overline {v}$

Die Ableitungsmatrizen enthalten die Koeffizienten der entsprechenden Finiten Differenzen. Damit lässt sich die Matrixschreibweise der Orr-Sommerfeld-Gleichung aufstellen:

$\overline {A}_{0}\cdot \overline {v}+\overline {A}_{2}\cdot \overline {v}_{{yy}}+\overline {A}_{4}\cdot \overline {v}_{{yyyy}}=\omega \left(\overline {B}_{0}\cdot \overline {v}+\overline {B}_{2}\cdot \overline {v}_{{yy}}\right)$

$\left[\overline {A}_{0}+\overline {A}_{2}\cdot \overline {D}_{2}+\overline {A}_{4}\cdot \overline {D}_{4}\right]\cdot \overline {v}=\omega \left[\overline {B}_{0}+\overline {B}_{2}\cdot \overline {D}_{2}\right]\cdot \overline {v}$

$\overline {A}\cdot \overline {v}=\omega \cdot \overline {B}\cdot \overline {v}$

Bei Matrizen ${\overline {A}}$ , $\overline {A}_{2}$ , $\overline {A}_{4}$ sind ausschließlich die Hauptdiagonalen mit den entsprechenden Termen der Orr-Sommerfeld-Gleichung belegt. So besteht etwa die Hauptdiagonale der Matrix $\overline {A}_{2}$ aus den Elementen $a_{2}(j,j)=\alpha \cdot U0(j)-i{\frac {2\alpha ^{2}}{Re}}$ , wobei j der y-Index ist. Das Eigenwertproblem der resultierenden Matrizen ${\overline {A}}$ und ${\overline {B}}$ kann dann mittels numerischer Bibliotheken, z.B. LAPACK gelöst werden.

Ist man an der räumlichen Lösung interessiert, so besteht die Möglichkeit, den betrachteten Eigenwert auf das räumliche Problem zu iterieren. Hierzu wird der Imaginärteil $\alpha _{i}$ so iteriert, dass der Imaginärteil $\omega_i$ dieses Eigenwerts zu null wird. Die dafür notwendige Bestimmung des Eigenwertes einer nur leicht abgeänderten Matrix kann z.B. mit der Wielandt-Iteration erfolgen. Die eigentliche Iteration kann dann z.B. mittels des Newton-Verfahrens oder ähnlicher Methoden erfolgen.

Der große Vorteil des Matrixverfahrens besteht darin, dass man ohne Vorgabe von Schätzwerten das gesamte Spektrum erhält. Allerdings ist es relativ rechenaufwändig, vor allem bei der Iteration auf die räumliche Lösung.

Schießverfahren

Beim Schießverfahren wird das ursprüngliche Randwertproblem als ein Anfangswertproblem betrachtet und die fehlenden Anfangsbedingungen durch eine Iteration bestimmt. Hierzu sind die ursprünglichen Differentialgleichungen in ein System von Differentialgleichungen erster Ordnung umzuwandeln. Für die zweidimensionale Orr-Sommerfeld-Gleichung ergibt sich beispielhaft folgendes Gleichungssystem:

$z_{0}={\hat {v}}$

${\frac {\partial z_{0}}{\partial y}}=z_{1}$

${\frac {\partial z_{1}}{\partial y}}=z_{2}$

${\frac {\partial z_{2}}{\partial y}}=z_{3}$

${\frac {\partial z_{3}}{\partial y}}=\left(-\omega \alpha ^{2}+\alpha ^{3}u_{0}+\alpha u_{{0yy}}-{\frac {i\alpha ^{4}}{Re}}\right)\cdot z_{0}+\left(\omega -\alpha u_{0}+{\frac {i2\alpha ^{2}}{Re}}\right)\cdot z_{2}$

Dieses Gleichungssystem ist nun über der -Richtung zu integrieren. Für den Freistromrand (Parallelströmung außerhalb der Grenzschicht) lassen sich folgende Fundamentallösungen der Form ${\hat {v}}(y)\sim e^{{\lambda \cdot y}}$ herleiten:

reibungsfreie Fundamentallösung:

$\lambda _{{1,2}}=\pm \alpha$

viskose Fundamentallösung:

$\lambda _{{3,4}}=\pm {\sqrt {\alpha ^{2}+i\cdot Re(\omega +\alpha u_{0})}}$

Aufgrund des Abklingens der Störungen nach außen existiert nur eine physikalisch sinnvolle Lösung pro Fundamentallösung. Für das Beispiel einer Grenzschicht mit der Wand bei y=0 ergibt sich, dass aufgrund der Bedingung

${\hat {v}}(y=\infty )={\hat {v}}_{{y}}(y=\infty )=0$

nur die negativen Fundamentallösungen $\lambda _{{2,4}}<0$ physikalisch sinnvoll sind. Durch die Wahl der Integrationsrichtung ausgehend vom Fernfeld zur Wand hin lässt sich erreichen, dass die unphysikalischen Lösungen $\lambda _{{1,3}}>0$ gedämpft sind.

Das grundsätzliche Vorgehen für die räumliche Lösung z.B. einer Grenzschicht sieht dann so aus, dass die Frequenz $\omega _{r}$ und im allgemeinen dreidimensionalen Fall die spannweitige Wellenzahl $\gamma$ vorgegeben werden. Ausgehend vom Freistromrand werden die beiden Fundamentallösungen $\lambda _{{2,4}}$ für eine geschätzte Wellenzahl $\alpha$ entlang der -Richtungen integriert. Dabei ist es in der Regel erforderlich, eine Orthonormalisierung durchzuführen, um ein Anwachsen der anderen Fundamentallösung zu unterdrücken.

Die Linearkombination der beiden Fundamentallösungen und die Wellenzahl $\alpha$ sind dann so zu iterieren, sodass die Wandrandbedingungen: ${\hat {v}}(y=0)={\hat {v}}_{{y}}(y=0)=0$ erfüllt werden. Bei der Integration über der -Richtung ist natürlich die Stabilitätsgrenze des verwendeten Verfahrens einzuhalten.

Der Hauptvorteil des Schießverfahrens ist seine Schnelligkeit, da dabei keine großen Gleichungssysteme zu lösen sind. Es eignet sich gleichermaßen für das zeitliche und das räumliche Problem. Nachteilig ist die Vorgabe von sinnvollen Schätzwerten, was insbesondere für hohe Machzahlen ein Problem darstellen kann. Die für die Grenzschicht dargestellte Vorgehensweise lässt sich nur bedingt auf andere Strömungen wie eine freie Scherschicht übertragen.

Literatur

H. Schlichting, K. Gersten: Grenzschicht-Theorie. 9. Auflage. Springer Verlag, Berlin 1997, ISBN 3-540-55744-X.
L. Mack: Boundary-layer linear stability theory. In: AGARD Special Course on Stability and Transition of Laminar Flow. AGARD R-709, 1984.

Siehe auch

Basierend auf einem Artikel in:

Wikipedia.de