BA Emilio Schmidt
Grundlagen
Gesucht ist ein Zusammenhang für die Kraft \(\vec{F}\), die auf eine von einem viskosen Fluid mit der Einstromgeschwindigkeit \(\vec{v}=-v\vec{e}_z\) umströmte Kugel mit dem Radius \(R\) wirkt.
Dabei werden hier stark laminare Strömungen betrachtet. Eine Größe, die die Laminarität bzw. Turbulenz einer Strömung charakterisiert ist die Reynolds-Zahl \(\text{Re}\). Sie ist eine dimensionslose Größe, welche definiert ist durch: \begin{align*} \text{Re}=\frac{\rho vL}{\mu} \end{align*}
Hierbei ist \(\rho\) die Dichte des Fluids, \(v\) die Geschwindigkeit des Fluids, mit der es auf das angeströmte Objekt zufließt, \(L\) die charakteristische Länge des angeströmten Körpers und \(\mu\) die dynamische Viskosität des Fluids. Stark laminare Strömungen werden dann durch \(\text{Re}\ll 1\) beschrieben.
Analytische Herleitung der Stokes'schen Formel
Die auf die Kugel wirkende Kraft lässt sich in zwei separate Beiträge unterteilen. Einerseits kann eine Kraft aufgrund eines Druckgradienten um die Kugel wirken, andererseits kann eine Kraft infolge von Scherungseffekten an der Kugeloberfläche auftreten.
Die Kraft auf die Kugel lässt sich dann über den Spannungstensor \(\underline{\underline{\sigma}}\) berechnen: \begin{align*} \vec{F}=\oint_{\partial V} \underline{\underline\sigma}\cdot\text d\vec A \end{align*}
Der Spannungstensor beschreibt dabei die Verteilung von Kräften bzw. Spannungen in einem Material. Er gibt an jedem Punkt eines Kontinuums die internen Kräfte bzw. Spannungen an, die aufgrund von äußeren Krafteinwirkungen wirken. Der Spannungstensor ist ein Tensor 2. Stufe, dessen Komponenten \(\sigma_{ij}\) die Spannungen, die aufgrund von Kräften in Richtung der \(j\)-ten Koordinatenachse auf eine Fläche mit einer Normalen in Richtung \(i\)-ten Koordinatenachse wirken, repräsentieren. Im Allgemeinen ist der Spannungstensor gegeben durch: \begin{align*} \underline{\underline{\sigma}}= -[p-\zeta(\vec{\nabla}\cdot\vec{u})]\mathbb{1} +\mu\big[\vec{\nabla}\vec{u}+(\vec{\nabla}\vec{u})^T -\frac{2}{3}(\vec{\nabla}\cdot\vec{u})\mathbb{1}\big] \end{align*}
Dabei ist \(p\) das Druckfeld, \(\zeta\) die Volumenviskosität, \(\vec{u}\) das Geschwindigkeitsfeld und \(\mathbb{1}\) der Einheitstensor. Um den Spannungstensor bestimmen zu können, müssen also zunächst das Geschwindigkeitsfeld \(\vec{u}\) und das Druckfeld \(p\) bestimmt werden. Dafür müssen die Kontinuitätsgleichung und die Navier-Stokes-Gleichung gelöst werden, welche im Allgemeinen gegeben sind durch: \begin{align*} &\frac{\partial\rho}{\partial t}+\vec{\nabla}\cdot(\rho\vec{u})=0 \\ &\rho\bigg(\frac{\partial\vec{u}}{\partial t} +(\vec{u}\cdot\vec{\nabla})\vec{u}\bigg)= -\vec{\nabla}p+\mu\Delta\vec{u} +\frac{\mu}{3}\vec{\nabla}(\vec{\nabla}\cdot\vec{u}) +\vec{f}_{\text{ext}} \end{align*}
Dabei ist \(\vec{f}_{\text{ext}}\) eine externe Kraftdichte. Diese Gleichungen sind jedoch im Allgemeinen nicht analytisch zu lösen, weshalb zusätzliche Annahmen über das System zu treffen sind. Dabei wird hier davon ausgegangen, dass es sich bei vorliegendem Fluid um ein Newtonsches Fluid handelt, woraus folgt, dass die dynamische Viskosität \(\mu=\text{const.}\) ist. Zusätzlich folgt dann nach der Stokes'schen Hypothese, dass die Volumenviskosität für Newtonsche Fluide vernachlässigt werden kann: \begin{align*} \zeta=0 \end{align*}
Außerdem lässt sich, aufgrund der stark laminaren Strömung annehmen, dass das Fluid inkompressibel ist, da die Dichteänderungen bei solch langsamen Strömungen sehr gering ist. Dabei bedeutet Inkompressibilität, dass die Dichte räumlich homogen und zeitlich konstant ist: \begin{align*} \rho=\text{const.} \end{align*}
Zudem kann, wegen der stark laminaren Strömung angenommen werden, dass diese stationär ist. Dies liegt daran, da die Reynolds-Zahl dem Verhältnis von Trägheitskräften, die versuchen das Fluid in Bewegung zu halten, zu viskosen Kräften, die Bewegungsunterschiede innerhalb des Fluides versuchen zu dämpfen, entspricht. Für Reynolds-Zahlen \(\text{Re}\ll 1\) dominieren die viskosen Kräfte, was dazu führt, dass das Fluid in geordneten Schichten, ohne Querbewegung fließt. Aus der Stationärität des Fluids folgt, dann dass jede Zeitableitung von Größen, die dieses Fluid beschreiben, verschwinden.
Darüber hinaus kann, wegen obiger Argumentation, der Konvektionsterm \(\rho(\vec{u}\cdot\vec{\nabla})\vec{u}\) vernachlässigt werden, da eben dieser die Trägheitskräfte beschreibt.
Mit der Annahme, dass keine externen Kräfte vorliegen, haben die zu lösenden Gleichungen dann folgende Form: \begin{align*} &\vec{\nabla}\cdot\vec{u}=0 \\ &\vec{\nabla}p=\mu\Delta\vec{u} \end{align*}
Mit der Kontinuitätsgleichung für inkompressible Newtonsche Fluide folgt dann für den Spannungstensor: \begin{align*} \underline{\underline{\sigma}}= -p\mathbb{1} +\mu\big[\vec{\nabla}\vec{u}+(\vec{\nabla}\vec{u})^T\big] \end{align*}
Die zu lösenden Gleichungen sind nun analytisch lösbar. Es wird sich als sinnvoll herausstellen zunächst das Geschwindigkeitsfeld zu bestimmen. Dafür ist es nützlich die Rotation der Navier-Stokes-Gleichung zu betrachten. Denn mit der Tatsache, dass Gradientenfelder rotationsfrei sind, folgt: \begin{align*} \vec{\nabla}\times(\Delta\vec{u})=\Delta(\vec{\nabla}\times\vec{u})=\vec{0} \end{align*}
Aufgrund dessen, dass für das Geschwindigkeitsfeld im Unendlichen \(\vec{u}(r\rightarrow\infty)\) gelten muss, lässt sich das Geschwindigkeitsfeld folgendermaßen schreiben: \begin{align*} \vec{u}=\vec{u}^\prime+\vec{v} \end{align*}
Dabei muss jedoch \(\vec{u}^\prime\) im Unendlichen verschwinden.
Einsetzen in die Kontinuitätsgleichung liefert dann mit \(\vec{v}=\text{const.}\): \begin{align*} \vec{\nabla}\cdot\vec{u}^\prime=0 \end{align*}
Dies entspricht der notwendigen Bedingung für die Existenz eines Vektorpotentials \(\vec{A}\) zu dem Vektorfeld \(\vec{u}^\prime\). Daher lässt sich \(\vec{u}^\prime\) folgendermaßen schreiben: \begin{align*} \vec{u}^\prime=\vec{\nabla}\times\vec{A} \end{align*}
Da das Geschwindigkeitsfeld \(\vec{u}\) und somit auch \(\vec{u}^\prime\) und die Rotation eines polaren Vektors ein axialer Vektor ist und umgekehrt, muss das Vektorpotential ein axialer Vektor sein. Außerdem darf das Vektorpotential nur vom Ortsvektor \(\vec{r}=r\vec{e}_r$\), sowie der Einstromgeschwindigkeit \(\vec{v}\) abhängen, da dieses über das Geschwindigkeitsfeld \(\vec{u}^\prime\) bestimmt wird und \(\vec{u}^\prime\) eben nur abhängig vom Ortsvektor und der Einstromgeschwindigkeit ist. Da sowohl der Ortsvektor, als auch die Einstromgeschwindigkeit polare Vektoren sind, muss das Vektorpotential in irgendeiner Form aus dem Vektorprodukt dieser beiden bestehen. Es bietet sich also an folgenden Ansatz für das Vektorpotential zu wählen: \begin{align*} \vec{A}=f^\prime(r)\vec{e}_r\times\vec{v}=\vec{\nabla}f(r)\times\vec{v} \end{align*}
Dabei ist \(f(r)\) eine skalare Funktion, welche nur vom Betrag des Ortsvektors \(r\) abhängt. Damit folgt für das Geschwindigkeitsfeld \(\vec{u}^\prime\): \begin{align*} &\phantom{\leftrightarrow\,,} \vec{u}^\prime=\vec{\nabla}\times(\vec{\nabla}f(r)\times\vec{v})\,, \;\text{mit}\;\vec{v}=\text{const.} \Rightarrow\vec{\nabla}\times(f\vec{v})=\vec{\nabla}f\times\vec{v}\\ &\Leftrightarrow \vec{u}^\prime=\vec{\nabla}\times\vec{\nabla}\times(f\vec{v}) \end{align*}
Einsetzen liefert für das Geschwindigkeitsfeld \(\vec{u}\): \begin{align*} \vec{u}=\vec{\nabla}\times\vec{\nabla}\times(f\vec{v})+\vec{v} \end{align*}
Einsetzen in die Navier-Stokes-Gleichung und ausnutzen der Relation \(\vec{\nabla}\times\vec{u}=-\Delta(\vec{\nabla}\times(f\vec{v}))\) liefert:\begin{align*} &\phantom{\leftrightarrow\,,} \Delta^2(\vec{\nabla}\times(f\vec{v}))=\vec{0}\,, \;\text{mit}\;\vec{v}=\text{const.} \Rightarrow\vec{\nabla}\times(f\vec{v})=\vec{\nabla}f\times\vec{v}\\ &\Leftrightarrow \Delta^2(\vec{\nabla}f\times\vec{v})=\vec{0}\\ &\Leftrightarrow \Delta^2(\vec{\nabla}f)\times\vec{v}=\vec{0} \end{align*}
Aufgrund dessen, dass obiger Ausdruck für alle Einstromgeschwindigkeiten gelten muss, muss insbesondere gelten: \begin{align*} \Delta^2(\vec{\nabla}f)=\vec{0} \end{align*}
Eine erste Integration liefert dann: \begin{align*} \Delta^2f=\text{const.} \end{align*}
Da das Geschwindigkeitsfeld \(\vec{u}^\prime\) im Unendlichen verschwinden muss, müssen auch dessen Ableitungen im Unendlichen verschwinden. Aufgrund dessen, dass \(\Delta^2f\) die vierten Ableitungen von \(f\) enthält und das Geschwindigkeitsfeld \(\vec{u}^\prime\) die zweiten Ableitungen von \(f\) enthält, muss die Konstante Null sein. Damit folgt: \begin{align*} \Delta^2f=0 \end{align*}
Da die Funktion \(f(r)\) ausschließlich vom Betrag des Ortsvektors \(r\) abhängt, reduziert sich der Laplace-Operator zu: \begin{align*} \Delta^2f=\frac{1}{r^2}\frac{\partial}{\partial r} \bigg(r^2\frac{\partial}{\partial r}\bigg)\Delta f=0 \end{align*}
Damit ist \(\Delta f\) allgemein gegeben durch: \begin{align*} \Delta f=\frac{2a}{r}+A \,, \quad a,A\in\mathbb{R} \end{align*}
Damit das Geschwindigkeitsfeld \(\vec{v}^\prime\) im Unendlichen verschwindet, müssen alle additiven Konstanten gleich Null gewählt werden. Einsetzen liefert: \begin{align*} \Delta f=\frac{2a}{r} \end{align*}
Analoges Vorgehen zu oben liefert dann für die Funktion \(f(r)\): \begin{align*} f(r)=ar+\frac{b}{r} \,, \quad a,b\in\mathbb{R} \end{align*}
Einsetzen in \(\vec{u}=\vec{u}=\vec{\nabla}\times\vec{\nabla}\times(f\vec{v})+\vec{v}\) liefert dann: \begin{align*} \vec{u}=\vec{v}-a\frac{\vec{v}+(\vec{v}\cdot\vec{e}_r)\vec{e}_r}{r} +b\frac{3(\vec{v}\cdot\vec{e}_r)\vec{e}_r-\vec{v}}{r^3} \end{align*}
Die Integrationskonstanten \(a\) und \(b\) lassen sich nun über die Randbedingung, dass das Geschwindigkeitsfeld auf der Kugeloberfläche verschwinden muss bestimmen. Es muss also gelten: \begin{align*} &\phantom{\leftrightarrow\,,} \vec{u}(r=R)=\vec{0}\\ &\Leftrightarrow \vec{0}=\vec{v}-a\frac{\vec{v}+(\vec{v}\cdot\vec{e}_r)\vec{e}_r}{r} +b\frac{3(\vec{v}\cdot\vec{e}_r)\vec{e}_r-\vec{v}}{r^3}\\ &\Leftrightarrow \vec{0}=\bigg(-\frac{a}{R}-\frac{b}{R^3}+1\bigg)\vec{v} +\bigg(-\frac{a}{R}+\frac{3b}{R^3}\bigg)(\vec{v}\cdot\vec{e}_r)\vec{e}_r \end{align*}
Aufgrund dessen, dass die obige Gleichung für alle \(\vec{v}\) und alle \(\vec{e}_r\) gelten muss, müssen deren Koeffizienten verschwinden. Daher sind die Integrationskonstanten gegeben durch: \begin{align*} a=\frac{3}{4}R\quad b=\frac{1}{4}R^3 \end{align*}
Damit ist das Geschwindigkeitsfeld endgültig gegeben durch: \begin{align*} \vec{u}=\vec{v}-\frac{3}{4}R\frac{\vec{v}+(\vec{v}\cdot\vec{e}_r)\vec{e}_r}{r} -\frac{1}{4}R^3\frac{\vec{v}-3(\vec{v}\cdot\vec{e}_r)\vec{e}_r}{r^3} \end{align*}
Die jeweiligen Komponenten sind dann gegeben durch: \begin{align*} & u_r=v\cos(\theta)\bigg(1-\frac{3R}{2r}+\frac{R^3}{2r^3}\bigg)\\ & u_\theta=v\sin(\theta)\bigg(1-\frac{3R}{4r}-\frac{R^3}{4r^3}\bigg)\\ & u_\varphi=0 \end{align*}
Um nun das Druckfeld berechnen zu können, wird der Ausdruck \(\vec{u}=\vec{\nabla}\times\vec{\nabla}\times(f\vec{v})+\vec{v}\) in die Navier-Stokes-Gleichung eingesetzt und die Relation \(\vec{\nabla}\times\vec{\nabla}\times(f\vec{v})=\vec{\nabla}(\vec{\nabla}f\cdot\vec{v})\vec{v}\Delta f\) ausgenutzt: \begin{align*} &\phantom{\leftrightarrow\,,} \vec{\nabla}p=\mu\Delta(\vec{\nabla}(\vec{\nabla}f\cdot\vec{v})-\vec{v}\Delta f) \,,\;\text{mit}\;\Delta^2f=0\\ &\Leftrightarrow \vec{\nabla}p=\mu\Delta(\vec{\nabla}(\vec{\nabla}f\cdot\vec{v}))\\ &\Leftrightarrow \vec{\nabla}p=\vec{\nabla}(\mu\Delta(\vec{\nabla}f\cdot\vec{v})) \end{align*}
Eine erste Integration liefert: \begin{align*} &\phantom{\leftrightarrow\,,} p=\mu\Delta(\vec{\nabla}f\cdot\vec{v})+p_0\\ &\Leftrightarrow p=\mu \vec{v}\cdot\vec{\nabla}(\Delta f)+p_0 \,,\;\text{mit}\; f=\frac{3}{4}Rr+\frac{1}{4}\frac{R^3}{r}\\ &\Leftrightarrow p=p_0-\frac{3}{2r^2}\mu R\vec{v}\cdot\vec{e}_r\\ &\Leftrightarrow p=p_0+\frac{3}{2r^2}\mu vR\cos(\theta) \end{align*}
Dabei ist \(p_0\) der Druck des Fluides im Unendlichen. Nun kann also die Kraft auf die Kugel berechnet werden. Dabei gilt für ein Oberflächenelement \(\text d\vec{A}\) einer Kugel mit Radius \(r=R\), in Kugelkoordinaten: \begin{align*} d\vec{A}=R^2\sin(\theta) d\varphi d\theta\vec{e}_r \end{align*}
Einsetzen in obige Gleichung für die Kraft liefert dann: \begin{align*} \vec{F}=\int_{0}^{\pi}\int_{0}^{2\pi} \underline{\underline\sigma}(r=R,\theta,\varphi)\cdot\vec{e}_r R^2\sin(\theta) d\varphi d\theta \end{align*}
Auswerten des Matrix-Vektor-Produktes liefert dann, dass nur die \(rr\)-, \(\theta r\)- und \(\varphi r\)-Komponenten des Spannungstensors multipliziert mit entsprechenden Basisvektoren übrig bleiben. Es gilt also: \begin{align*} \vec{F}=\int_{0}^{\pi}\int_{0}^{2\pi} [\sigma_{rr}\vec{e}_r+\sigma_{\theta r}\vec{e}_{\theta} +\sigma_{\varphi r}\vec{e}_{\varphi}]\Big|_{r=R} R^2\sin(\theta) d\varphi d\theta \end{align*}
Die \(rr\)-Komponente des Spannungstensors ist dann gegeben durch: \begin{align*} \sigma_{rr}=-p+2\mu(\vec{\nabla}\vec{u})_{rr} \end{align*}
\(\vec{\nabla}\vec{u}_{rr}\) ist hier die \(rr\)-Komponente des Vektorgradienten \(\vec{\nabla}\vec{u}\), welche gegeben ist durch: \begin{align*} &\phantom{\leftrightarrow\,,} \vec{\nabla}\vec{u}_{rr}=\frac{\partial u_r}{\partial r} \,,\;\text{mit}\; u_r=v\cos(\theta)\bigg(1-\frac{3R}{2r}+\frac{R^3}{2r^3}\bigg)\\ &\Leftrightarrow \vec{\nabla}\vec{u}_{rr}=0 \end{align*}
Mit \(p=p_0+\frac{3}{2r^2}\mu vR\cos(\theta)\) folgt dann für \(rr\)-Komponente des Spannungstensors auf der Kugeloberfläche: \begin{align*} \sigma_{rr}\Big|_{r=R}=-p_0-\frac{3}{2R}\mu v\cos(\theta) \end{align*}
Für die \(\theta r\)-Komponente gilt dann gemäß obiger Gleichung: \begin{align*} \sigma_{\theta r}=\mu\bigg[ \vec{\nabla}\vec{u}_{\theta r} +\vec{\nabla}\vec{u}_{r\theta}\bigg] \end{align*}
Einsetzen der jeweiligen Komponenten des Vektorgradienten und auswerten an der Kugeloberfläche liefert: \begin{align*} &\phantom{\leftrightarrow\,,} \sigma_{\theta r}\Big|_{r=R}= \mu\bigg[ \frac{\partial u_{\theta}}{\partial r}\Big|_{r=R} +\frac{1}{R} \frac{\partial u_r}{\partial\theta}\Big|_{r=R} -\frac{u_{\theta}}{R}\Big|_{r=R}\bigg]\\ &\Leftrightarrow \sigma_{\theta r}\Big|_{r=R}=\frac{3}{2R}\mu v\sin(\theta) \end{align*}
Für die \(\varphi r\)-Komponente gilt gemäß obiger Gleichung: \begin{align*} \sigma_{\varphi r}=\mu\bigg[ \vec{\nabla}\vec{u}_{\varphi r} +\vec{\nabla}\vec{u}_{r\varphi}\bigg] \end{align*}
Einsetzen der jeweiligen Komponenten des Vektorgradienten und auswerten an der Kugeloberfläche liefert: \begin{align*} &\phantom{\leftrightarrow\,,} \sigma_{\varphi r}\Big|_{r=R}= \bigg[ \frac{\partial u_{\varphi}}{\partial r}\Big|_{r=R} +\frac{1}{R\sin(\theta)} \frac{\partial u_r}{\partial\varphi}\Big|_{r=R} -\frac{u_\varphi}{R}\Big|_{r=R}\bigg]\\ &\Leftrightarrow \sigma_{\varphi r}\Big|_{r=R}=0 \end{align*}
Damit folgt für die Kraft auf die umströmte Kugel: \begin{align*} &\phantom{\leftrightarrow\,,} \vec{F}=\int_{0}^{\pi}\int_{0}^{2\pi} [-p_0\vec{e}_{r}-\frac{3}{2R}\mu v\cos(\theta)\vec{e}_{r} +\frac{3}{2R}\mu v\sin(\theta)\vec{e}_{\theta}]\Big|_{r=R} R^2\sin(\theta) d\varphi d\theta\\ &\Leftrightarrow \vec{F}=6\pi\mu R\vec{v} \end{align*}
Oseen'sche Näherung
Bei der Herleitung der Stokes'schen Formel wurden Annahmen getroffen, die für sehr kleine Reynolds-Zahlen \(\text{Re}\ll 1\) gültig sind. Oseen schaffte es nun eine erste Näherung für kleine Reynolds-Zahlen \(\text{Re}<1\) zu formulieren, indem er den Konvektionsterm, in linearer Form, berücksichtigt. Aufgrund dessen, dass in der Oseen'schen Näherung immer noch vergleichsweise kleine Reynolds-Zahlen betrachtet werden, kann weiterhin von einer stationären Strömung ausgegangen werden. Des weiteren kann, wegen der vergleichsweise kleinen Reynolds-Zahlen weiterhin von einem inkompressiblen Fluid ausgegangen werden. Mit diesen Annahmen folgt für die Navier-Stokes-Gleichung: \begin{align*} \rho(\vec{u}\cdot\vec{\nabla})\vec{u}= -\vec{\nabla}p+\mu\Delta\vec{u} \end{align*}
Um nun den Konvektionsterm linearisieren zu können, wird analog zu oben davon ausgegangen, dass sich das Geschwindigkeitsfeld folgendermaßen schreiben lässt: \begin{align*} \vec{u}=\vec{u}^\prime+\vec{v}\,, \;\text{mit}\;\vec{u}\prime(r\rightarrow\infty)=\vec{0} \end{align*}
Einsetzen in den Konvektionsterm liefert: \begin{align*} &\phantom{\leftrightarrow\,,} (\vec{u}\cdot\vec{\nabla})\vec{u} =((\vec{u}^\prime+\vec{v})\cdot\vec{\nabla}) (\vec{u}^\prime+\vec{v})\\ &\Leftrightarrow (\vec{u}\cdot\vec{\nabla})\vec{u} =(\vec{u}^\prime\cdot\vec{\nabla})\vec{u}^\prime +(\vec{u}^\prime\cdot\vec{\nabla})\vec{v} +(\vec{v}\cdot\vec{\nabla})\vec{u}^\prime +(\vec{v}\cdot\vec{\nabla})\vec{v}\,, \;\text{mit}\;\vec{v}=\text{const.}\\ &\Leftrightarrow (\vec{u}\cdot\vec{\nabla})\vec{u} =(\vec{u}^\prime\cdot\vec{\nabla})\vec{u}^\prime +(\vec{v}\cdot\vec{\nabla})\vec{u}^\prime \end{align*}
Da hier weiterhin kleine Reynolds-Zahlen betrachtet werden, kann davon ausgegangen werden, dass die Störung \(\vec{u}^\prime\) klein ist, sodass Terme zweiter Ordnung, wie \((\vec{u}^\prime\cdot\vec{\nabla})\vec{u}^\prime\) vernachlässigt werden können. Damit folgt für den linearisierten Konvektionsterm: \begin{align*} (\vec{u}\cdot\vec{\nabla})\vec{u} =(\vec{v}\cdot\vec{\nabla})\vec{u}^\prime \end{align*}
Damit folgt für die Navier-Stokes-Gleichung: \begin{align*} \rho(\vec{v}\cdot\vec{\nabla})\vec{u}^\prime= -\vec{\nabla}p+\mu\Delta\vec{u}^\prime \end{align*}
Die zu lösenden Gleichungen sind also gegeben durch: \begin{align*} &\vec{\nabla}\cdot\vec{u}^\prime=0 \\ &\rho(\vec{v}\cdot\vec{\nabla})\vec{u}^\prime= -\vec{\nabla}p+\mu\Delta\vec{u}^\prime \end{align*}
Lösen der Gleichungen und berechnen der Kraft gemäß \(\vec{F}=\oint_{\partial V}\underline{\underline{\sigma}}\cdot d\vec{A}\) liefert dann eine Korrektur erster Ordnung für die Stokes'sche Formel: \begin{align*} \vec{F}=6\pi\mu R\vec{u}\bigg(1+\frac{3}{8}\mathrm{Re}\bigg) \end{align*}
Numerische Berechnung der Kraft auf die Kugel
\(\newcommand{\et}{ {\large\mathbb 1} }\)
Im Unterschied zu der Annahme, die für die analytische Herleitung getroffen wurde, simuliert das verwendete Simulationsprogramm PLUTO kompressible Newtonsche Fluide. Damit verschwindet die Divergenz des Geschwindigkeitsfeldes nicht mehr, weshalb der Spannungstensor für solche Fluide gegeben ist durch: \begin{align*} \underline{\underline\sigma}= p\mathbb{1} +\mu\bigg[\vec{\nabla}\vec{u}+(\vec{\nabla}\vec{u})^T -\frac{2}{3}(\vec{\nabla}\cdot\vec{u})\mathbb{1}\bigg] \end{align*}
Für die Kraft sind dabei jedoch nur die \(rr\)-, \(\theta r\)- und \(\phi r\)-Komponenten des Spannungstensors von Interesse. Die \(rr\)-Komponente ist dabei gegeben durch: \begin{align*} \sigma_{rr}=-p+\mu\bigg[2\vec{\nabla}\vec{u}_{rr} -\frac{2}{3}\vec{\nabla}\cdot\vec{u}\bigg] \end{align*}
\(\vec{\nabla}\vec{u}_{rr}\) ist hier die \(rr\)-Komponente des Vektorgradienten \(\vec{\nabla}\vec{u}\), welche gegeben ist durch: \begin{align*} \vec{\nabla}\vec{u}_{rr}=\frac{\partial u_r}{\partial r} \end{align*}
Für die Divergenz des Geschwindigkeitsfeldes gilt allgemein in Kugelkoordinaten: \begin{align*} \vec{\nabla}\cdot\vec{u}= \frac{1}{r^2}\frac{\partial r^2 u_r}{\partial r} +\frac{1}{r\sin(\theta)} \frac{\partial\sin(\theta)u_{\theta}}{\partial\theta} +\frac{1}{r\sin(\theta)} \frac{\partial u_{\varphi}}{\partial\varphi} \end{align*}
Aufgrund dessen, dass das Fluid homogen einströmt und eine Kugel umströmt, hat das Geschwindigkeitsfeld um die Kugel Zylindersymmetrie. Daher fällt die azimutale Ableitung bei Betrachtung des Geschwindigkeitsfeldes um die Kugel weg. Außerdem folgt aus der vorgegebenen Anfangsbedingung, dass das Fluid entlang der z-Achse einströmt, zusammen mit der Kugelsymmetrie der Kugel, dass das Geschwindigkeitsfeld und auch sämtliche anderen relevanten Vektorfelder, die dieses Problem beschreiben, keine azimutale Komponente aufweisen können. Damit folgt für die Divergenz des Geschwindigkeitsfeldes:
\begin{align*} \vec{\nabla}\cdot\vec{u}= \frac{1}{r^2} \frac{\partial r^2 u_r}{\partial r} +\frac{1}{r\sin(\theta)} \frac{\partial\sin(\theta)u_{\theta}}{\partial\theta} \end{align*}
Einsetzen in die \(rr\)-Komponente und auswerten an der Kugeloberfläche liefert: \begin{align*} \sigma_{rr}\Big|_{r=R}=-p +\mu\bigg[ 2\frac{\partial u_r}{\partial r}\Big|_{r=R} -\frac{2}{3}\frac{1}{R^2} \frac{\partial r^2 u_r}{\partial r}\Big|_{r=R} +\frac{1}{R\sin(\theta)} \frac{\partial\sin(\theta)u_{\theta}}{\partial\theta}\Big|_{r=R}\bigg] \end{align*}
Die Ableitungen lassen sich zusätzlich mittels der Produktregel vereinfachen. Es folgt: \begin{align*} \sigma_{rr}\Big|_{r=R}=-p +\mu\bigg[ \frac{4}{3}\frac{\partial u_r}{\partial r}\Big|_{r=R} -\frac{4}{3R}u_r\Big|_{r=R} -\frac{2}{3R}\cot(\theta)u_{\theta}\Big|_{r=R} -\frac{2}{3R}\frac{\partial u_{\theta}} {\partial\theta}\Big|_{r=R}\bigg] \end{align*}
Mit den Randbedingungen folgt, dass für das Geschwindigkeitsfeld gelten muss: \begin{align*} u_r\Big|_{r=R}=u_{\theta}\Big|_{r=R}=0 \end{align*}
Außerdem impliziert das Verschwinden der Geschwindigkeitsvektoren auf der Kugeloberfläche und dessen Kugelsymmetrie, dass zusätzlich die polaren Ableitungen auf der Oberfläche verschwinden: \begin{align*} \frac{\partial u_r}{\partial\theta}\Big|_{r=R}= \frac{\partial u_{\theta}}{\partial\theta}\Big|_{r=R}=0 \end{align*}
Damit folgt für die \(rr\)-Komponente des Spannungstensors, ausgewertet an der Kugeloberfläche: \begin{align*} \sigma_{rr}\Big|_{r=R}=-p+\frac{4}{3}\mu \frac{\partial u_r}{\partial r}\Big|_{r=R} \end{align*}
Für die \(\theta r\)-Komponente gilt analog zu oben: \begin{align*} \sigma_{\theta r}=\mu\bigg[ \vec{\nabla}\vec{u}_{\theta r} +\vec{\nabla}\vec{u}_{r\theta}\bigg] \end{align*}
Einsetzen der jeweiligen Komponenten des Vektorgradienten und auswerten an der Kugeloberfläche liefert: \begin{align*} \sigma_{\theta r}\Big|_{r=R}= \mu\bigg[ \frac{\partial u_{\theta}}{\partial r}\Big|_{r=R} +\frac{1}{R} \frac{\partial u_r}{\partial\theta}\Big|_{r=R} -\frac{u_{\theta}}{R}\Big|_{r=R}\bigg] \end{align*}
Ausnutzen der Randbedingungen liefert dann: \begin{align*} \sigma_{\theta r}\Big|_{r=R}= \mu\frac{\partial u_{\theta}}{\partial r}\Big|_{r=R} \end{align*}
Für die \(\varphi r\)-Komponente gilt gemäß obiger Gleichung: \begin{align*} \sigma_{\varphi r}=\mu\bigg[ \vec{\nabla}\vec{u}_{\varphi r} +\vec{\nabla}\vec{u}_{r\varphi}\bigg] \end{align*}
Die jeweiligen Komponenten des Vektorgradienten sind im Allgemeinen gegeben durch: \begin{align*} &\vec{\nabla}\vec{u}_{\varphi r}= \frac{\partial u_{\varphi}}{\partial r}\\ &\vec{\nabla}\vec{u}_{r\varphi}= \frac{1}{r\sin(\theta)} \frac{\partial u_r}{\partial\varphi} -\frac{u_\varphi}{r} \end{align*}
Analog zu oben fallen hier auch sowohl die azimutalen Komponenten als auch die azimutalen Ableitungen weg, weshalb für die \(\varphi r\)-Komponente des Spannungstensors folgt: \begin{align*} \sigma_{\varphi r}=0 \end{align*}
Einsetzen in das Integral zur Kraftberechnung und auswerten des Integrals über \(\varphi\) liefert: \begin{align*} &\phantom{\leftrightarrow,,} \vec{F}=2\pi R^2\int_{0}^{\pi}\bigg[\bigg( -p +\frac{4}{3}\mu\frac{\partial u_r} {\partial r}\Big|_{r=R}\bigg)\vec{e}_r +\mu \frac{\partial u_{\theta}}{\partial r}\Big|_{r=R} \vec{e}_{\theta}\bigg]\sin(\theta)\text d\theta\\ &\Leftrightarrow \vec{F}=2\pi R^2 \bigg[ -\int_{0}^{\pi}p\vec{e}_r\sin(\theta)d\theta +\frac{4}{3}\int_{0}^{\pi}\mu \frac{\partial u_r}{\partial r}\Big|_{r=R} \vec{e}_r\sin(\theta)\text d\theta +\int_{0}^{\pi}\mu \frac{\partial u_{\theta}}{\partial r}\Big|_{r=R} \vec{e}_{\theta}\sin(\theta)\text d\theta\bigg] \end{align*}
Da im vorliegenden Problem insbesondere die Kraft auf die Kugel in \(z\)-Richtung von Interesse ist, wird nun die Projektion des Kraftvektors auf die z-Achse bestimmt: \begin{align*} F_z=\vec{e}_z\cdot\vec{F} \end{align*}
Aufgrund dessen, dass die \(z\)-Achse eine feste Richtung hat und der Linearität der Integrationsoperation, kann der Basisvektor \(\vec{e}_z\) in das Integral hineingezogen werden. Damit folgt für \(F_z\): \begin{align*} &F_z=2\pi R^2 \bigg[ -\int_{0}^{\pi}p\vec{e}_r\cdot\vec{e}_z \sin(\theta)d\theta +\frac{4}{3}\int_{0}^{\pi}\mu \frac{\partial u_r}{\partial r}\Big|_{r=R} \vec{e}_r\cdot\vec{e}_z\sin(\theta)\text d\theta\\ &\phantom{\leftrightarrow,,} +\int_{0}^{\pi}\mu \frac{\partial u_{\theta}}{\partial r}\Big|_{r=R} \vec{e}_{\Theta}\cdot\vec{e}_z \sin(\theta)\text d\theta\bigg] \end{align*}
Für die Berechnung der jeweiligen Skalarprodukte ist es praktisch den Basisvektor \(\vec{e}_z\) in der hier genutzten sphärischen Basis darzustellen: \begin{align*} \vec{e}_z=\cos(\theta)\vec{e}_r -\sin(\theta)\vec{e}_{\theta} \end{align*}
Damit folgt für \(F_z\): \begin{align*} &F_z=2\pi R^2 \bigg[ -\int_{0}^{\pi}p\cos(\theta)\sin(\theta)d\theta +\frac{4}{3}\int_{0}^{\pi}\mu \frac{\partial u_r}{\partial r}\Big|_{r=R} \cos(\theta)\sin(\theta)\text d\theta\\ &\phantom{\leftrightarrow,,} -\int_{0}^{\pi}\mu \frac{\partial u_{\theta}}{\partial r}\Big|_{r=R} \sin^2(\theta)\text d\theta\bigg] \end{align*}
Die Kraft in \(x\)- bzw. \(y\)-Richtung lässt sich analog berechnen. Da jedoch die Basisvektoren \(\vec{e}_x\) bzw. \(\vec{e}_y\) in sphärischer Darstellung eine \(\varphi\)-Abhängigkeit aufweisen, kann die Integration über \(\varphi\) nicht so leicht durchgeführt werden, wie es oben geschehen ist. Für die Kraft in \(x\)-Richtung muss dann folgendes Integral ausgewertet werden: \begin{align*} &F_x=R^2 \bigg[ -\int_{0}^{\pi}\int_{0}^{2\pi} p(R,\theta)\sin^2(\theta)\cos(\varphi)d\varphi d\theta +\frac{4}{3}\int_{0}^{\pi}\int_{0}^{2\pi}\mu(R,\theta) \frac{\partial u_r}{\partial r}\Big|_{r=R} \sin^2(\theta)\cos(\varphi)d\varphi d\theta\\ &\phantom{\leftrightarrow,,} +\int_{0}^{\pi}\int_{0}^{2\pi}\mu(R,\theta) \frac{\partial u_{\theta}}{\partial r}\Big|_{r=R} \sin(\theta)\cos(\theta)\cos(\varphi)d\varphi d\theta\bigg] \end{align*}
Auswerten der Integrale über \(\varphi\) liefert jedoch, wie zu erwarten war, für \(F_x\): \begin{align*} F_x=0 \end{align*}
Analoges Vorgehen liefert für \(F_y\): \begin{align*} F_y=0 \end{align*}
Numerische Integration
Da PLUTO die Daten nur für die jeweiligen Zellzentren ausgibt, werden insbesondere keine Daten an der Kugeloberfläche ausgegeben. Diese Daten sind jedoch gemäß der obigen Formel zur Berechnung der Kraft erforderlich. Das Geschwindigkeitsfeld an der Kugeloberfläche ist jedoch bereits aus den Randbedingungen bekannt: \begin{align*} \vec{u}(r=R,\theta)=\vec{0}\quad \forall\;\theta\in[0,\pi] \end{align*}
Der Druck auf der Kugeloberfläche ist jedoch unbekannt und muss daher extrapoliert werden. Eine lineare Extrapolation bietet sich an, da sie einen Fehler zweiter Ordnung aufweist, welcher mit den Fehlern der verwendeten Methoden zur Bestimmung der Ableitungen der Felder übereinstimmt. Zudem ist eine lineare Extrapolation ausreichend, da bei einer ausreichend hohen Gitterauflösung die Zellzentren der Zellen in der ersten Zellschale um die Kugeloberfläche nicht weit von dieser entfernt sind. Für die lineare Extrapolation des Drucks auf die Kugeloberfläche $p(r=R)=p_0$ werden die beiden nächsten Datenpunkte \((r_1,p_1)\) und \((r_2,p_2)\) benötigt. Der extrapolierte Druck ist dann gegeben durch: \begin{align*} p_0=p_2+\frac{R-r_2}{r_1-r_2}(p_1-p_2) \end{align*}
Aufgrund dessen, dass die Daten jedoch keine kontinuierlichen Größen sind, kann das für die Kraftberechnung benötigte Integral nur näherungsweise berechnet werden. Obiges Integral ist dabei ein Integral über den Polarwinkel \(\theta\) über eine Funktion \(f(\theta)\), welche nur von diesem Winkel abhängt. Um das Integral also zu nähern, wird der obige Integrand \(f(\theta_i)\) in jeder Zelle für den Winkel $\theta_i$ ausgewertet und dann Gewichtet mit der jeweiligen Breite der Zelle \(\Delta\theta_i\) in polarer Richtung, über alle \(N_\theta\) aufsummiert: \begin{align*} \int_{0}^{\pi}f(\theta)d\theta \rightarrow\sum_{i}f(\theta_i)\Delta\theta_i \end{align*}
Aufgrund dessen, dass die polare Aufteilung der Zellen in vorliegendem Gitter uniform gewählt wurde gilt: \begin{align*} \Delta\theta_i=\Delta\theta_j=\Delta\theta\quad \forall\;i,j\in[0,N_\theta-1] \end{align*}
Damit folgt: \begin{align*} \int_{0}^{\pi}f(\theta)d\theta \rightarrow\sum_{i}f(\theta_i)\Delta\theta \end{align*}
Kompressibilität des Fluides
Wie oben angemerkt, simuliert das verwendete Simulationsprogramm PLUTO kompressible Fluide. Die Stokes'sche Formel wurde jedoch unter der Annahme, dass es sich bei dem vorliegenden Fluid um ein inkompressibles Fluid handelt, hergeleitet. Um die Vergleichbarkeit mit der analytisch hergeleiteten Formel sicherzustellen, muss die Kompressibilität des Fluids \(\beta\) während jeder durchgeführten Simulation konstant gehalten werden. In der Software PLUTO kann die Kompressibilität jedoch nicht manuell konstant gesetzt werden, da sie eine abgeleitete Größe aus anderen Parametern ist. Die Kompressibilität ist dabei per Definition gegeben durch: \begin{align*} \beta=\frac{1}{\rho}\frac{\partial\rho}{\partial p} \end{align*}
Damit lässt sich die Kompressibilität als Funktion der Schallgeschwindigkeit \(c_S\) in dem Fluid ausdrücken. Denn die Schallgeschwindigkeit ist im Allgemeinen gegeben durch: \begin{align*} c_S=\sqrt{\frac{\partial p}{\partial\rho}} \end{align*}
Damit folgt für die Kompressibilität als Funktion der Schallgeschwindigkeit: \begin{align*} \beta=\frac{c_S^2}{\rho} \end{align*}
Aufgrund dessen, dass bei stark laminaren Strömungen keine großen Dichteänderungen auftreten, kann die Dichte \(\rho\) des Fluides folgendermaßen geschrieben werden: \begin{align*} \rho=\rho_\infty+\rho^\prime \end{align*}
Dabei ist \(\rho_\infty\) die Dichte des Fluides im Unendlichen und \(\rho^\prime\) die entsprechende Störung, für die gelten muss: \begin{align*} \rho^\prime\ll 1 \end{align*}
Damit folgt für die Kompressibilität: \begin{align*} \beta=\frac{c_S^2}{\rho_\infty+\rho^\prime} \approx\frac{c_S^2}{\rho_\infty} \end{align*}
Damit also die Kompressibilität des Fluides konstant bleibt, müssen also die Schallgeschwindigkeit des Fluides und die Dichte des Fluides im Unendlichen konstant gehalten werden: \begin{align*} \beta=\text{const.}\Leftrightarrow c_S=\text{const.}\land\rho_\infty=\text{const.} \end{align*}
Natürliche Einheiten
In der Physik und insbesondere bei der Durchführung von Experimenten ist es üblich, numerische Ergebnisse in SI-Einheiten auszudrücken. SI-Einheiten bieten eine universelle Grundlage, die es ermöglicht, experimentelle Daten und theoretische Ergebnisse konsistent und vergleichbar zu machen. Allerdings kann die Verwendung von SI-Einheiten bei theoretischen Berechnungen zu komplexen und schwer überschaubaren Gleichungen führen. Um diese Berechnungen zu vereinfachen und vergleichbarer zu gestalten, werden sogenannte natürliche Einheiten verwendet. Natürliche Einheiten sind speziell gewählte Maßeinheiten, die an die charakteristischen Größen der untersuchten physikalischen Systeme angepasst sind. Ein wesentlicher Vorteil dieser Einheiten ist die Vereinfachung physikalischer Gleichungen. Durch die geeignete Wahl von natürlichen Einheiten können viele physikalische Konstanten auf \(1\) gesetzt werden, was die Anzahl der in den Gleichungen auftretenden Konstanten erheblich reduziert. Das Ziel bei der Verwendung natürlicher Einheiten besteht darin, die betreffenden Gleichungen dimensionslos zu machen. Für das gegebene Problem werden die relevanten Größen in Einheiten der Dichte im Unendlichen \(\rho_\infty\), der dynamischen Viskosität \(\mu\) und des Objektradius $R$ gemessen. Bei der Wahl der natürlichen Einheiten ist es entscheidend, dass diese die physikalisch relevanten Dimensionen korrekt reproduzieren. Die Massen-, Längen- und Zeiteinheiten in natürlichen Einheiten sind dann gegeben durch: \begin{align*} \rho_\infty R^3, \; R, \; \frac{\rho_\infty R^2}{\mu} \end{align*}
Mit der Definition der Reynolds-Zahl ergibt sich für die Einstromgeschwindigkeit in den gewählten natürlichen Einheiten: \begin{align*} \frac{v}{\mu/(\rho_\infty R)} = \frac{\text{Re} \mu / (2 \rho_\infty R)}{\mu / (\rho_\infty R)} = \frac{\text{Re}}{2} \end{align*}
Da die Reynolds-Zahl eine dimensionslose Größe ist, ist die Einstromgeschwindigkeit in den gewählten natürlichen Einheiten dimensionslos. Damit ergibt sich die Stokes'sche Formel in den gewählten natürlichen Einheiten zu: \begin{align*} \frac{F}{\mu^2 / \rho_\infty} = 3 \pi \text{Re} \end{align*}
Diese ist nun, wie gefordert, eine dimensionslose Größe. \colorbox{yellow}{Sache mit Schallgeschwindigkeit}
Natürliche Einheiten (Lothar)
Aus den vier System-Parametern \(v,~R,~\rho\) und \(\mu\) lassen sich zwei sinnvolle Einheitensätze bilden:
"zähe" Einheiten
Damit ist der Satz \(R,~\rho\) und \(\mu\) gemeint, Basiseinheiten sind $$ \rho R^3~,\quad R~,\quad\frac{R^2}{\rho\mu} $$ für Masse, Länge und Zeit; und \(\mu/(\rho R)\) bzw. \(\mu^2/\rho\) die Geschwindigkeits- bzw. Krafteinheit. Die Reynoldszahl stellt dann gemäß $$ \frac{v}{\mu/(\rho R)} = \frac{\text{Re}}{2} $$ die Geschwindigkeit ein, und im Stokes-Regime gilt $$ \frac{F}{\mu^2/\rho} = 3\pi\text{Re}~. $$
"träge" Einheiten
Damit ist der Satz \(R,~v\) und \(\rho\) gemeint, Basiseinheiten sind dann $$ \rho R^3~,\quad R~,\quad\frac{R}{v} $$ für Masse, Länge und Zeit; und \(\rho vR\) bzw. \(\rho v^2R^2\) die Viskositäts- bzw. Krafteinheit. Die Reynoldszahl stellt dann gemäß $$ \frac{\mu}{\rho vR} = \frac{2}{\text{Re}} $$ die Viskosität ein, und im Stokes-Regime gilt $$ \frac{F}{\rho v^2R^2} = 6\pi\frac{2}{\text{Re}}~. $$
NB: Die Zahlenwerte der einheitenbildenden Parameter sind 1 und sollten auch im Code so gewählt werden.