BA Pia Lüdecke

Erste Tests

Setup von Emilio zeigt eine mit regelmäßigen Abständen in Stufen abnehmende \(E_\text{kin}\). Warum? (Mit besserer Gitterauflösung keine Stufen)

\(E_\text{kin}\) steigt nicht wie befürchtet exponentiell, sondern erreicht ein etwas höher liegendes Plateau. Möglicherweise erst dort stationärer Zustand?

Der stationäre Zustand scheint ungefähr bei \(t_{end}=1e6\) erreicht zu werden.

Grundlagen

Näherungen

In einem Paper aus dem Jahr 2016 fand A. Morrison einen Fit, der den Widerstandskoeffizienten für die Strömung um eine Kugel bis hin zu \(Re<10^6\) approximiert:

\[C_D \approx \frac{24}{Re} + \frac{2,6 \cdot (\frac{Re}{5})} {1+(\frac{Re}{5})^{1,52}} + \frac{0,411 \cdot (\frac{Re}{2,63 \cdot 10^5})^{-7,94}}{1+(\frac{Re}{2,63 \cdot 10^5})^{-8}} + \frac{0,25 \cdot (\frac{Re}{10^6})}{1+(\frac{Re}{10^6})}\]

Für laminare Strömung und bis zu \(Re < 2 \cdot 10^5\) kann nach Kaskas eine kürzere Version einer Näherungslösung genutzt werden:

\[C_D = \frac{24}{Re} + \frac{4}{\sqrt{Re}} +0,4 \]

Im Bereich sehr kleiner Reynoldszahlen \(Re<1\) können die letzten beiden Terme vernachlässigt werden und es folgt das Gesetz von Stokes:

\[C_D = \frac{24}{Re}\]

Numerische Berechnung der Kraft auf die Kugel

Da zur Berechnung der Kraft auf eine Kugel nur die \(rr-\),\(\theta r-\) und \(\varphi r-\)Komponente des Spannungstensors übrig bleiben, hat die zu berechnende Kraft die Form:

\[\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 \]

Der Spannungstensor wird um die Volumenviskosität \(\zeta\) erweitert. Damit hat er die Form:

\[\sigma = -\Big[p- \zeta \vec{\nabla} \cdot \vec{u}\Big] \mathbb{1} + \mu \Big[ \vec{\nabla} \vec{u} +(\vec{\nabla} \vec{u})^T - \frac{2}{3} (\vec{\nabla} \cdot \vec{u}) \mathbb{1} \Big] \]

Zunächst soll die \(rr\)-Komponente betrachtet werden, die gegeben ist durch:

\[\sigma_{rr} = -\Big[p- \zeta \vec{\nabla} \cdot \vec{u}\Big] + \mu\Big[ 2\vec{\nabla} \vec{u}_{rr} - \frac{2}{3} (\vec{\nabla} \cdot \vec{u}) \Big] \]

Mit der Divergenz in Kugelkoordinaten, die gegeben ist durch:

\[\vec{\nabla} \cdot \vec{u} = \frac{1}{r^2} \partial_r (r^2 u_r) + \frac{1}{r \sin{(\theta})} \partial_{\theta}(\sin{(\theta)} u_{\theta}) + \frac{1}{r \sin{(\theta})} \partial_{\phi} u_{\phi} \]

ergibt sich daher:

\[\sigma_{rr} = -\Big[p- \zeta ( \frac{2}{r} u_r + \partial_r u_r+ \frac{1}{r} \partial_{\theta} u_{\theta} + \frac{\cot{(\theta)}}{r} u_{\theta} )\Big] + \mu \Big[\frac{4}{3} \partial_r u_r - \frac{4}{3r} u_r - \frac{2}{3r} u_{\theta} \cot{(\theta)} - \frac{2}{3r} \partial_{\theta} u_{\theta} \Big] \]

Dabei wurde bereits berücksichtigt, dass bei dem betrachteten Problem aufgrund der Zylindersymmetrie die azimutale Ableitung entfällt. Unter Berücksichtigung der Randbedingungen gilt weiterhin für das Geschwindigkeitsfeld:

\[u_r \big \vert _{r=R} = u_{\theta} \big \vert _{r=R} = 0\]

Aufgrund der no-slip-condition muss weiterhin gelten, dass die polaren Ableitungen auf der Kugelberfläche verschwinden:

\[\partial_{\theta} u_r \big \vert _{r=R} = \partial_{\theta} u_{\theta} \big \vert _{r=R} = 0\]

Damit folgt das finale Ergebnis für die \(rr\)-Komponente:

\[\sigma_{rr} \big \vert _{r=R} = -p + \zeta \partial_r u_r \big \vert _{r=R} + \mu\frac{4}{3} \partial_r u_r \big \vert _{r=R} \]

Analog werden nun auch die \(\theta r\)- und \(\varphi r\)-Komponente berechnet. Da die Volumenviskosität dort nicht eingeht, ergibt sich das bekannte Ergebnis:

\[\sigma_{\theta r} \big \vert _{r=R} = \mu \partial_r u_{\theta} \big \vert _{r=R} \] \[\sigma_{\varphi r} \big \vert _{r=R} = 0 \]

Eingesetzt in die oben gegebene Formel folgt somit für die Kraft auf die Kugel:

\[\vec{F}=2\pi R^2\int_{0}^{\pi}\Big[(-p + \zeta \partial_r u_r \big \vert _{r=R} + \mu\frac{4}{3} \partial_r u_r \big \vert _{r=R})\vec{e}_r+ \mu \partial_r u_{\theta} \big \vert _{r=R}\vec{e}_{\theta} \Big] \sin(\theta) d\theta \]

Da erwartungsgemäß die \(y\)- und \(x\)-Komponente der Kraft auf die Kugel null werden, wird im Folgenden nur die \(z\)-Komponente der Kraft genauer betrachtet:

\[\vec{F_z}=2\pi R^2\Big[(-\int_{0}^{\pi} p \vec{e}_r \cdot \vec{e}_z \sin(\theta) d\theta + \int_{0}^{\pi} \zeta \partial_r u_r \big \vert _{r=R} \vec{e}_r \cdot \vec{e}_z \sin(\theta) d\theta + \frac{4}{3} \int_{0}^{\pi} \mu\ \partial_r u_r \big \vert _{r=R} \vec{e}_r \cdot \vec{e}_z \sin(\theta) d\theta + \int_{0}^{\pi} \mu \partial_r u_{\theta} \big \vert _{r=R}\vec{e}_{\theta} \cdot \vec{e}_z \sin(\theta) d\theta \Big] \]

Um weiter mit der Rechnung fortzufahren ist es notwendig \(\vec{e_z}\) in der sphärischen Basis darzustellen:

\[ \vec{e_z} = \cos{(\theta)} \vec{e_r} - \sin{(\theta)} \vec{e_{\theta}} \]

Damit ergibt sich das finale Ergebnis:

\[\vec{F_z}=2\pi R^2\Big[(-\int_{0}^{\pi} p \cos{(\theta)} \sin{(\theta)} d\theta + \int_{0}^{\pi} \zeta \partial_r u_r \big \vert _{r=R} \cos{(\theta)} \sin(\theta) d\theta + \frac{4}{3} \int_{0}^{\pi} \mu\ \partial_r u_r \big \vert _{r=R} \cos{(\theta)} \sin(\theta) d\theta - \int_{0}^{\pi} \mu \partial_r u_{\theta} \big \vert _{r=R} \sin(\theta)^2 d\theta \Big] \]