Turbulenzmodell

Ein Turbulenzmodell dient in der numerischen Strömungssimulation zur Schließung des zugrundeliegenden Gleichungssystems.

Da Turbulenz sowohl räumlich als auch zeitlich auf sehr unterschiedlichen und vor allem auch sehr kleinen Skalen stattfindet, werden zur korrekten Auflösung aller Phänomene extrem feine Gitter und Zeitschritte benötigt, wie dies in der direkten numerischen Simulation (DNS) geschieht.

Die DNS ist aufgrund ihres extremen Rechenaufwands auf absehbare Zeit auf niedrige (und damit für praktische Anwendungen häufig wenig interessante) Reynolds-Zahlen beschränkt. Aus diesem Grund werden unterschiedliche Strategien zur Verringerung des Rechenaufwandes verwendet. Je mehr Annahmen dabei getroffen werden, desto geringer wird in der Regel die Rechenzeit und desto größer werden die Unsicherheiten bezüglich des Ergebnisses.

Statistische Modellierung

Die derzeit am weitesten verbreitete Modellierungsstrategie ist die statistische Modellierung (auch: Reynolds Averaged Navier Stokes (RANS)). Dabei wird die turbulente Strömung im einfachsten Fall als zeitlicher Mittelwert ${\overline {\bullet }}$ und Varianz ${\bullet '}$ der Geschwindigkeit $u_{i}$ und des Drucks $p$ modelliert.

Führt man diese Modellannahme in die Navier-Stokes-Gleichungen ein und mittelt diese zeitlich, so entstehen die Reynoldsgemittelten Navier-Stokes-Gleichungen. Diese enthalten mit dem Reynolds'schen Spannungstensor

${-\tau _{ij}}=\rho \,{\overline {u_{i}'u_{j}'}}=\rho {\begin{pmatrix}{\overline {u_{1}'u_{1}'}}&{\overline {u_{1}'u_{2}'}}&{\overline {u_{1}'u_{3}'}}\\{\overline {u_{2}'u_{1}'}}&{\overline {u_{2}'u_{2}'}}&{\overline {u_{2}'u_{3}'}}\\{\overline {u_{3}'u_{1}'}}&{\overline {u_{3}'u_{2}'}}&{\overline {u_{3}'u_{3}'}}\end{pmatrix}}$

einen zusätzlichen Term mit zusätzlichen Variablen $\rho \,{\overline {u_{i}'u_{j}'}}$, den Reynoldsspannungen. Die Diagonalenelemente des Tensors stellen Normalspannungen dar, während die restlichen Elemente Schubspannungen sind. Das Gleichungssystem ist damit nicht mehr geschlossen. Die Schließung gelingt durch zusätzliche Annahmen für die Komponenten des Reynolds'schen Spannungstensors in Form von Gleichungen. Diese zusätzlichen Gleichungen werden Turbulenzmodell genannt. Da Turbulenz noch im Wesentlichen unverstanden ist, beruhen diese meist auf Heuristiken. Zur Validierung dienen Daten aus Experimenten. Man unterscheidet zwischen Null-, Ein- und Zweigleichungsmodellen, sowie Schließungsansätzen 2. Ordnung.

Wirbelviskositätsmodelle

Bei den Wirbelviskositätsmodellen wird der Reynolds'sche Spannungstensor nach der Näherung von Boussinesq approximiert. Hierbei werden die Reynoldsspannungen in Analogie zu den durch molekulare Viskosität hervorgerufenen Spannungen behandelt:

$\rho {\overline {u_{i}'u_{j}'}}=-\mu _{t}\left({\frac {\partial {\bar {u}}_{i}}{\partial x_{j}}}+{\frac {\partial {\bar {u}}_{j}}{\partial x_{i}}}\right)+{\frac {2}{3}}\rho k\delta _{ij}$.

Die Größe $\mu _{t}=\rho \cdot \nu _{t}$ wird turbulente Wirbelviskosität genannt und beschreibt die Erhöhung der Viskosität durch turbulente Schwankungsbewegungen. In der Regel überschreitet $\mu _{t}$ die molekulare Viskosität deutlich. Division durch die Dichte ergibt die turbulente kinematische Wirbelviskosität, $\nu _{t}$. Die Wurzel aus der turbulenten kinetischen Energie $k={\overline {u_{i}'u_{i}'}}/2$ stellt ein typisches Geschwindigkeitsmaß der turbulenten Schwankungsbewegungen dar. Das Symbol $\delta _{ij}$ bezeichnet das Kronecker-Delta. Der Term ${\frac {2}{3}}\rho k\delta _{ij}$ ist ein „turbulenter Druckterm“, der notwendig ist, um die Gleichung auch für Normalspannungen anwenden zu können.

Aus Dimensionsgründen lässt sich die turbulente kinematische Wirbelviskosität $\nu _{t}$ mit einem turbulenten Längenmaß $L_{t}$ und einem turbulenten Geschwindigkeitsmaß $U_{t}$ gemäß $\nu _{t}\sim L_{t}\cdot U_{t}$ ausdrücken. Durch den Boussinesq-Ansatz gelingt das Schließen der Reynoldsgemittelten Navier-Stokes Gleichungen durch die Bestimmung der Wirbelviskosität $\nu _{t}$ bzw. der damit assoziierten turbulenten Längen- und Geschwindigkeitsmaße.

Die Wirbelviskositätsmodelle werden nach der Anzahl der unabhängigen Turbulenzvariablen, die zur Berechnung von $L_{t}$ und $U_{t}$ benutzt werden, unterschieden.

Nullgleichungsmodelle

Algebraische oder Nullgleichungsmodelle verwenden zur Schließung nur algebraische Beziehungen. Hierzu zählen das Baldwin-Lomax-Modell und das Turbulenzmodell nach Cebeci und Smith.

Eingleichungsmodelle

Eingleichungsmodelle benutzen eine zusätzliche Transportgleichung zur Bestimmung von $\nu _{t}$. Das verbreitetste Eingleichungsmodell stammt von Spalart und Allmaras, das eine zusätzliche Transportgleichung für die an die turbulente Viskosität angelehnte Hilfsgröße ${\tilde {\nu }}$ einführt. Außer in Wandnähe stimmt ${\tilde {\nu }}$ mit der turbulenten Viskosität $\nu _{t}$ überein:

$\rho {\frac {\partial {\tilde {\nu }}}{\partial t}}+\rho {\frac {\partial \left({\tilde {\nu }}u_{i}\right)}{\partial x_{i}}}={\frac {1}{\sigma _{\tilde {\nu }}}}\left[{\frac {\partial }{\partial x_{j}}}\left\{\left(\mu +\rho {\tilde {\nu }}\right){\frac {\partial {\tilde {\nu }}}{\partial x_{j}}}\right\}+C_{b2}\rho \left({\frac {\partial {\tilde {\nu }}}{\partial x_{j}}}\right)^{2}\right]-C_{\omega 1}\rho f_{\omega }\left({\frac {\tilde {\nu }}{d}}\right)^{2}+C_{b1}\rho {\tilde {S}}{\tilde {\nu }}$

Die beiden Terme hinter der eckigen Klammer beschreiben die Turbulenzdestruktion und die Turbulenzproduktion. Nachteilig ist bei diesem Turbulenzmodell die Unfähigkeit, schnelle Änderungen im turbulenten Längenmaß, wie sie beim Übergang einer Grenzschicht in eine freie Scherschicht auftreten, richtig vorherzusagen.

Zweigleichungsmodelle

Zweigleichungs-Turbulenzmodelle sind ein Schließungsansatz, der aus der Lösung zweier gekoppelter Transportgleichungen besteht. Man unterscheidet die Modelle anhand der verwendeten Turbulenzgrößen. Zwei große Gruppen sind z. B. die $k{\text{-}}\varepsilon$-Turbulenzmodelle und die $k{\text{-}}\omega$-Turbulenzmodelle.

Standard k-ε-Turbulenzmodell

Das $k{\text{-}}\varepsilon$-Turbulenzmodell ist ein weitverbreitetes Zweigleichungsmodell. Es beschreibt mit zwei partiellen Differentialgleichungen die Entwicklung der turbulenten kinetischen Energie $k$ und der isotropen Dissipationsrate $\varepsilon =\nu {\overline {(\partial u_{i}'/\partial x_{k})(\partial u_{i}'/\partial x_{k})}}$. Die Gleichungen lauten:

$\rho {\frac {\partial k}{\partial t}}+\rho {\bar {u}}_{j}{\frac {\partial k}{\partial x_{j}}}=C_{\mu }\rho \mu _{t}\left({\frac {\partial {\bar {u}}_{i}}{\partial x_{j}}}+{\frac {\partial {\bar {u}}_{j}}{\partial x_{i}}}\right){\frac {\partial {\bar {u}}_{i}}{\partial x_{j}}}-\rho \varepsilon +{\frac {\partial }{\partial x_{j}}}\left[\left(\mu +{\frac {\mu _{t}}{\sigma _{k}}}\right){\frac {\partial k}{\partial x_{j}}}\right]$

und

$\rho {\frac {\partial \varepsilon }{\partial t}}+\rho {\bar {u}}_{j}{\frac {\partial \varepsilon }{\partial x_{j}}}=C_{\varepsilon 1}{\frac {\varepsilon }{k}}\tau _{ij}{\frac {\partial {\bar {u}}_{i}}{\partial x_{j}}}-C_{\varepsilon 2}{\frac {\varepsilon ^{2}}{k}}C_{\mu }\rho \mu _{t}\left({\frac {\partial {\bar {u}}_{i}}{\partial x_{j}}}+{\frac {\partial {\bar {u}}_{j}}{\partial x_{i}}}\right){\frac {\partial {\bar {u}}_{i}}{\partial x_{j}}}-C_{\varepsilon 2}\rho {\frac {\varepsilon ^{2}}{k}}+{\frac {\partial }{\partial x_{j}}}\left[\left(\mu +{\frac {\mu _{t}}{\sigma _{\varepsilon }}}\right){\frac {\partial \varepsilon }{\partial x_{j}}}\right]$

In die oben genannten Gleichungen sind einige, zum Teil erheblich vereinfachende Modellannahmen eingeflossen. Dies schränkt den Gültigkeitsbereich und somit den Anwendungsbereich deutlich ein. In den Gleichungen tauchen noch unbekannte Koeffizienten auf. Diese werden durch die Betrachtung einfacher Strömungsfelder ermittelt. Der Parameter $C_{\varepsilon 1}$ wird durch eine homogene Scherung im Gleichgewichtszustand kalibriert. Die Größe $C_{\varepsilon 2}$ folgt aus dem Abklingverhalten homogener Gitterturbulenz. Die turbulente Prandtlzahl $\sigma _{\varepsilon }$ ergibt sich aus einer Analyse des logarithmischen Bereiches einer ebenen turbulenten Wandgrenzschicht. Der Anisotropieparameter $C_{\mu }$ ergibt sich aus einer Dimensionsanalyse der Wirbelviskosität: $\nu _{t}\sim k\cdot (k/\varepsilon )$. Daraus folgt unmittelbar $\nu _{t}=C_{\mu }(k^{2}/\varepsilon )$. Die Betrachtung einer turbulenten Wandgrenzschicht liefert dann einen Wert für $C_{\mu }$.

Für das Standard $k{\text{-}}\varepsilon$-Modell findet man in der Literatur häufig:

${\begin{matrix}C_{\mu }&=&0{,}09\\C_{\varepsilon 1}&=&1{,}44\\C_{\varepsilon 2}&=&1{,}92\\\sigma _{\varepsilon }&=&1{,}3\\\sigma _{k}&=&1\\\end{matrix}}$

Die Art der Bestimmung der Konstanten bezeichnet die Strömungsfelder, in denen das Modell gute Übereinstimmung mit Messungen liefern sollte.

Nichtlineare k-ε-Turbulenzmodelle

Das Standard $k{\text{-}}\varepsilon$-Modell hat einige gravierende Nachteile. Die Normalspannungen werden durch die Boussinesq-Approximation des Reynolds'schen Spannungstensors in allen Raumrichtungen gleich groß berechnet. Dies darf nicht mit der klassischen Definition der Isotropie eines Tensors 2. Ordnung (Reynolds-Spannungstensor) verwechselt werden. Isotropie würde die Reynolds'schen Schubspannungen verschwinden lassen (isotrope Turbulenz), was jedoch nicht automatisch durch die Boussinesq-Approximation, welche Bestandteil des k-ε-Turbulenzmodells ist, abgebildet wird. Das bedeutet aber, dass Strömungsfelder, in denen der Geschwindigkeitsvektor in großem Maße von den Normalspannungen beeinflusst wird, nur ungenau abgebildet werden können. Dies ist der Fall in Ablösegebieten, Rezirkulationsbereichen und Sekundärströmungen. Ein Ausweg bietet eine Erweiterung der Boussinesq-Approximation. Dies führt zusätzliche nichtlineare Terme in die Modellgleichungen ein, die nichtlinear im Gradienten der mittleren Geschwindigkeit sind. Diese nichtlinearen Terme erlauben eine genauere Berechnung der Normalspannungen.

V2F-Turbulenzmodell

Die Turbulenz in der Nähe von Wänden ist gekennzeichnet durch Inhomogenität und Anisotropie. Die Zweigleichungsmodelle, wie $k{\text{-}}\varepsilon$ und $k{\text{-}}\omega$, benutzen in Wandnähe die Annahme homogener, isotroper Turbulenz. Dämpfungsfunktionen werden in diese Modelle eingefügt, um diese falschen Annahmen zu korrigieren. Dämpfungsfunktionen sind so gestaltet, dass bestimmte Lösungen vom Modell wiedergegeben werden können. In anderen Fällen werden falsche Vorhersagen getroffen. Das V2F-Turbulenzmodell ist eine Erweiterung des $k{\text{-}}\varepsilon$-Turbulenzmodells. Zusätzlich zu den Transportgleichungen für die turbulente kinetische Energie und die Dissipationsrate werden eine Gleichung für das Geschwindigkeitsmaß normal zur Wand, ${\overline {v'^{2}}}$ und deren mit $k$ normalisierte Produktionsrate $f$ gelöst. Die Gleichungen für $k$ und $\varepsilon$ sind identisch mit denen des Standard $k{\text{-}}\varepsilon$-Modells.

$\rho {\frac {\partial k}{\partial t}}+\rho {\bar {u}}_{j}{\frac {\partial k}{\partial x_{j}}}=C_{\mu }\rho \mu _{t}\left({\frac {\partial {\bar {u}}_{i}}{\partial x_{j}}}+{\frac {\partial {\bar {u}}_{j}}{\partial x_{i}}}\right)-\rho \varepsilon +{\frac {\partial }{\partial x_{j}}}\left[\left(\mu +{\frac {\mu _{t}}{\sigma _{k}}}\right){\frac {\partial k}{\partial x_{j}}}\right]$

und

$\rho {\frac {\partial \varepsilon }{\partial t}}+\rho {\bar {u}}_{j}{\frac {\partial \varepsilon }{\partial x_{j}}}=C_{\varepsilon 1}{\frac {\varepsilon }{k}}\tau _{ij}{\frac {\partial {\bar {u}}_{i}}{\partial x_{j}}}-C_{\varepsilon 2}{\frac {\varepsilon ^{2}}{k}}C_{\mu }\rho \mu _{t}\left({\frac {\partial {\bar {u}}_{i}}{\partial x_{j}}}+{\frac {\partial {\bar {u}}_{j}}{\partial x_{i}}}\right){\frac {\partial {\bar {u}}_{i}}{\partial x_{j}}}-C_{\varepsilon 2}\rho {\frac {\varepsilon ^{2}}{k}}+{\frac {\partial }{\partial x_{j}}}\left[\left(\mu +{\frac {\mu _{t}}{\sigma _{\varepsilon }}}\right){\frac {\partial \varepsilon }{\partial x_{j}}}\right]$

Für das wandnormale Geschwindigkeitsmaß wird die zusätzliche Gleichung

${\frac {\partial {\overline {v'^{2}}}}{\partial t}}+u\cdot \nabla {\overline {v'^{2}}}=kf_{22}-{\overline {v'^{2}}}{\frac {\varepsilon }{k}}+\nabla \left[\left(\nu +{\frac {\nu _{t}}{\sigma _{k}}}\right)\nabla {\overline {v'^{2}}}\right]$

formuliert. Der Term $kf_{22}$ stellt die Quelle für ${\overline {v'^{2}}}$ dar und kann als Umverteilung von Turbulenzintensität aus der strömungsparallelen Komponente interpretiert werden. Die nichtlokalen Effekte werden mathematisch durch eine elliptische Relaxationsgleichung für $f_{22}$ repräsentiert:

$L^{2}\nabla ^{2}f_{22}-f_{22}=\left(1-C_{1}\right){\frac {\left[{\frac {2}{3}}-{\frac {\overline {v'^{2}}}{k}}\right]}{T}}-{\frac {C_{2}}{k}}\nu _{t}\left({\frac {\partial u_{j}}{\partial x_{i}}}+{\frac {\partial u_{i}}{\partial x_{j}}}\right){\frac {\partial u_{j}}{\partial x_{i}}}$

Die im Model auftretenden Längen- und Zeitmaße sind:

$L=C_{L}l$

mit

$l^{2}=\max \left[{\frac {k^{3}}{\varepsilon ^{2}}},c_{\eta }^{2}\left({\frac {\nu ^{3}}{\varepsilon }}\right)^{\frac {1}{2}}\right]$

und

$T=\max \left[{\frac {k}{\varepsilon }},6\left({\frac {\nu }{\varepsilon }}\right)^{\frac {1}{2}}\right]$

Der Koeffizient $6$ im Ausdruck für $T$ wurde mit Hilfe direkter numerischer Simulation bestimmt. Die Wirbelviskosität ist gegeben durch:

$\nu _{t}=C_{\mu }{\overline {v'^{2}}}T$.

Die Modellkonstante $C_{\varepsilon 1}$ sollte je nach Wandabstand nach der Literatur zwischen $1{,}3$, weit entfernt von der Wand, und $1{,}55$, in einer anliegenden Grenzschicht, liegen. $C_{\varepsilon 1}$ wird mit der Gleichung

$C_{\varepsilon 1}=1{,}3+\left[{\frac {0{,}25}{1+\left({\frac {d}{2l}}\right)^{8}}}\right]$

zwischen diesen beiden Werten interpoliert. Die anderen Modellkonstanten sind gegeben mit:

${\begin{matrix}C_{\mu }&=&0{,}19\\C_{\varepsilon 2}&=&1{,}9\\C_{1}&=&1{,}4\\C_{2}&=&0{,}3\\C_{L}&=&0{,}3\\C_{\eta }&=&70{,}0\\\sigma _{\varepsilon }&=&1{,}3\\\sigma _{k}&=&1{,}0\\\end{matrix}}$

k-ω-Turbulenzmodell

Ein weiteres weitverbreitetes Zweigleichungs-Turbulenzmodell ist das von Wilcox angegebene $k{\text{-}}\omega$-Modell. Es werden hier eine Transportgleichung für $k$ und eine Transportgleichung für die charakteristische Frequenz, $\omega ={\frac {1}{C_{\mu }}}{\frac {\varepsilon }{k}}$, der energiedissipierenden Wirbel gelöst. Nach Wilcox lautet die Transportgleichung für $k$:

$\rho {\frac {\partial k}{\partial t}}+\rho {\bar {u}}_{j}{\frac {\partial k}{\partial x_{j}}}=C_{\mu }\rho \mu _{t}\left({\frac {\partial {\bar {u}}_{i}}{\partial x_{j}}}+{\frac {\partial {\bar {u}}_{j}}{\partial x_{i}}}\right){\frac {\partial {\bar {u}}_{i}}{\partial x_{j}}}-\beta ^{*}\rho k\omega +{\frac {\partial }{\partial x_{j}}}\left[\left(\mu +\sigma ^{*}\mu _{t}\right){\frac {\partial k}{\partial x_{j}}}\right]$

$\rho {\frac {\partial \omega }{\partial t}}+\rho {\bar {u}}_{j}{\frac {\partial \omega }{\partial x_{j}}}=\alpha {\frac {\omega }{k}}C_{\mu }\rho \mu _{t}\left({\frac {\partial {\bar {u}}_{i}}{\partial x_{j}}}+{\frac {\partial {\bar {u}}_{j}}{\partial x_{i}}}\right){\frac {\partial {\bar {u}}_{i}}{\partial x_{j}}}-\beta \rho \omega ^{2}+{\frac {\partial }{\partial x_{j}}}\left[\left(\mu +\sigma \mu _{t}\right){\frac {\partial \omega }{\partial x_{j}}}\right]$

$\beta ^{*}$ entspricht dem $C_{\mu }$ der $k{\text{-}}\varepsilon$-Modelle. Die Konstanten zur Schließung des Systems wurden in analoger Weise zum $k{\text{-}}\varepsilon$-Modell bestimmt und sind nach Wilcox gegeben mit:

$\alpha ={\frac {5}{9}},\beta ={\frac {3}{40}},\beta ^{*}={\frac {9}{100}},\sigma ={\frac {1}{2}},\sigma ^{*}={\frac {1}{2}}$

Das $k{\text{-}}\omega$-Modell reduziert das turbulente Längenmaß $L={\frac {k^{\frac {1}{2}}}{\omega }}$ in Wandnähe automatisch. Ein weiterer Vorteil besteht in der robusten Formulierung der viskosen Unterschicht. Nachteilig ist die Abhängigkeit des berechneten Grenzschichtrandes von der Freiströmbedingung für $\omega$, die vom Benutzer vorgegeben wird. Dieses Verhalten wird in der Literatur als „free stream“-Sensitivität bezeichnet.

k-ω-SST Turbulenzmodell

Das $k{\text{-}}\omega$-Modell bietet Vorteile in wandnahen Bereichen des Strömungsfeldes, wohingegen das $k{\text{-}}\varepsilon$-Modell in wandfernen Gebieten gute Resultate liefert. Die Vereinigung der Vorteile dieser beiden Modelle liefert das von Menter entwickelte SST-Turbulenzmodell.

Treten in der Strömung zusätzliche Phänomene (Verbrennung, Partikel, Tropfen, Überschall usw.) auf, so müssen auch die damit verbundenen Größen (bspw. Dichte, Temperatur, Massenbrüche etc.) gemittelt werden. In den dazugehörigen Transportgleichungen treten dabei analoge Schließungsprobleme auf.

Large Eddy Simulation

Anstelle der zeitlichen Mittelung wird bei der Large Eddy Simulation eine zeitliche und räumliche Tiefpassfilterung angewendet. Dies hat zur Folge, dass die großskaligen Phänomene transient simuliert werden, während der Beitrag der kleinskaligen Phänomene weiterhin modelliert werden.

Obwohl verwandte Modellierungsprobleme auftreten, verspricht die LES bei höherem Rechenaufwand eine bessere Beschreibung der Turbulenz als die statistischen Methoden, weil zumindest ein Teil der turbulenten Schwankungen wiedergegeben wird.

Detached Eddy Simulation

Die Detached Eddy Simulation (DES) wurde erstmals 1997 von P. Spalart veröffentlicht. Sie basiert in ihrer ursprünglichen Form auf dem Turbulenzmodell von Spalart-Allmaras (eine Transportgleichung), es wird aber auch an der Anwendung in Verbindung mit anderen Modellen geforscht.

Die DES ersetzt den Wandabstand, der als Variable im Spalart-Allmaras-Modell vorkommt, in wandfernen Bereichen durch die größte Weite einer Gitterzelle. Durch diese Formulierung lässt sich in den wandfernen Bereichen ein LES-ähnliches Verhalten der Rechnung erreichen. De facto erhält man so also eine RANS-Formulierung in der Grenzschicht und eine LES Formulierung in der freien Strömung, also das im jeweiligen Bereich am besten geeignete Verfahren (bezüglich Genauigkeit und Rechenaufwand).

Da RANS und LES unterschiedliche Anforderungen an das Gitter stellen, hat das Erstellen eines geeigneten, in entsprechende Zonen unterteilten Gitters einen großen Einfluss auf den Erfolg der Rechnung. Dasselbe gilt für die verwendeten numerischen Methoden. Diese sind aber meist gezwungenermaßen im gesamten Rechengebiet dieselben, was teilweise zu Kompromissen bezüglich der Genauigkeit führt.

Literatur

• Michael Breuer: Direkte Numerische Simulation und Large-Eddy Simulation turbulenter Strömungen auf Hochleistungsrechnern. 1. Auflage. Shaker, Aachen 2002, ISBN 3-8265-9958-6. 
• Herbert Oertel jr., Eckart Laurien: Numerische Strömungsmechanik. 2. Auflage. Vieweg, Braunschweig / Wiesbaden 2003, ISBN 3-528-03936-1. 
• Jochen Fröhlich: Large Eddy Simulation turbulenter Strömungen. 1. Auflage. Teubner, Wiesbaden 2006, ISBN 3-8351-0104-8. 
• Rüdiger Schwarze: CFD-Modellierung: Grundlagen und Anwendungen bei Strömungsprozessen, 1. Auflage. Springer, Berlin / Heidelberg, 2013, ISBN 978-3-642-24377-6.