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. Update: Bei \(t=0,2 \cdot 1e7\) macht die kinetische Energie einen weiteren Sprung. Dieser ist nur in der fünften Nachkommastelle bemerkbar, jedoch zeigt sich mit visit, dass sich an diesem Punkt die Geschwindigkeiten vertauschen (Geschwindigkeit um Kugel fällt plötzlich auf fast null, während die Geschwindigkeit im Kegel der Strömungsauftelung steigt). Im Folgenden soll geklärt werden, ob das Problem durch eine zu geringe Auflösung verursacht wird.

• Es ließ sich für \(Re=200\) kein stationärer Zustand finden, da unerwartete Turbulenzen auftraten. Im Folgenden wird die höchste Reynolds-Zahl auf \(Re=100\) reduziert.

• Mit \(Re=100\) stellt sich sehr schnell ein stationärer Zustand ein ( \(t_{end} \approx 500\) ). Dieser ist jedoch zunächst noch sehr abhängig von der gewählten Auflösung. Es werden nun unterschiedliche Auflösungen durchgetestet, um eine finale Konvergenz zu verifizieren. Update: Ab Auflösungen von n=9 zeigte sich, dass der stationäre Zustand bei \(t_{end} = 10^5\) doch noch nicht erreicht ist.

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 [1]:

\[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 des Weiteren für das Geschwindigkeitsfeld:

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

Aufgrund der no-slip-condition muss außerdem 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} \big(\zeta \partial_r u_r \big \vert _{r=R} -p + \frac{4}{3} \mu\ \partial_r u_r \big \vert _{r=R} \big)\cos{(\theta)} \sin(\theta) d\theta - \int_{0}^{\pi} \mu \partial_r u_{\theta} \big \vert _{r=R} \sin(\theta)^2 d\theta \Big] \]

Natürliche Einheiten

Da im Folgenden in allen Auswertungen der Widerstandskoeffizient anstelle der Widerstandskraft betrachtet werden soll, muss dieser noch in natürliche Einheiten umgerechnet werden. Im Allgemeinen gilt für den Widerstandskoeffizienten:

\[ C_D = \frac{2 F_D}{\rho v^2 A} \]

Wobei \(A\) die Stirnfläche ist und somit für eine umströmte Kugel \(A=\pi R^2\) gilt. Die Massen-, Längen- und Zeiteinheiten werden weiterhin in natürlichen Einheiten angegeben als:

\[ \rho_{\infty} R^3, R, \frac{R}{c_S} \]

Sind also gegeben durch den Kugelradius, die Dichte im Unendlichen und die Schallgeschwindigkeit. Im Falle einer symmetrischen Strömung (für Re=200 noch der Fall?) gilt für die Einströmgeschwindigkeit in natürlichen Einheiten:

\[ v^* = \frac{v}{c_S}= Ma\]

Die Bezugsfläche ist in natürlichen Einheiten gegeben als:

\[ A^* = \frac{A}{R^2} = \pi \]

Für die allgemeine Widerstandkraft in natürlichen Einheiten gilt:

\[ F_D^* = \frac{F_D}{R^2\rho_{\infty} c_S^2} \]

Eingesetzt in die Formel des Widerstandskoeffizienten folgt also:

\[ C_D = \frac{2 \cdot F_D^* \cdot (R^2\rho_{\infty} c_S^2)}{\rho A^* R^2 \cdot (v^* c_S)^2} \]

Vereinfachen und Kürzen des Ausdrucks liefert:

\[ C_D = \frac{2 F_D^* }{ v^{*2} \pi } \]