BA Emilio Schmidt

↑ up ↑

Stokes-Fluss um eine Kugel

Kraftberechnung

\(\newcommand{\et}{ {\large\mathbb 1} }\) \begin{gather} \vec F=\iint\limits_{\partial K} \underline{\underline\sigma}\cdot\text d\vec S\\ \text{und}\quad F_z = \vec e_z\cdot\iint\limits_{\partial K}\underline{\underline\sigma}\cdot\text d\vec S=\iint\limits_{\partial K}\vec e_z\cdot\underline{\underline\sigma}\cdot\text d\vec S \stackrel{\text{Kugel}}{=}\iint\limits_{\partial K}\vec e_z\cdot\underline{\underline\sigma}\cdot\vec e_r\,\text dS \end{gather} Bei Zylindersymmetrie und Verwendung von KK wird daraus \begin{align} F_z &= \int_0^\pi\int_0^{2\pi}(\vec e_r\cos\theta-\vec e_\theta\sin\theta)\cdot\underline{\underline\sigma}(r{=}R,\theta)\cdot\vec e_r\,R^2\sin\theta\,\text d\phi\,\text d\theta\\ &=2\pi R^2\int_0^\pi \Big(\sigma_{rr}(r{=}R,\theta)\cos\theta\sin\theta-\sigma_{\theta r}(r{=}R,\theta)\sin^2\theta\Big)\,\text d\theta \end{align}

Dabei lautet der Spannungstensor (s. auch Wikipedia) \begin{equation} \underline{\underline\sigma} = -\et(p - \zeta\vec\nabla\cdot\vec u) + \mu\Big(\vec\nabla\vec u + (\vec\nabla\vec u)^\text{T}-\et\frac{2}{3}\nabla\cdot\vec u\Big) \end{equation} und bei der Formulierung von \(\vec\nabla\vec u\) und \(\vec\nabla\cdot\vec u\) in KK hilft ebenfalls die Wikipedia. Die Komponenten des Einheitstensors \(\et\) schließlich bilden in jeder Basis die Einheitsmatrix.

Kraftberechnung (Emilio)

Die Kraft, welche auf die umströmte Kugel wirkt, ist gegeben durch folgendes Oberflächenintegral: \begin{align*} \vec{F}=\oint_{\partial V} \underline{\underline\sigma}\cdot\text d\vec A \end{align*}

Dabei gilt für ein Oberflächenelement \(\text d\vec{A}\) einer Kugel mit Radius \(r=R\), in Kugelkoordinaten: \begin{align*} \text d\vec{A}=R^2\sin(\theta)\text d\varphi\text d\theta\vec{e}_r \end{align*}

Einsetzen und weitere Berechnung in Kugelkoordinaten liefert: \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)\text d\varphi\text d\theta \end{align*}

Auswerten des Matrix-Vektor-Produktes liefert, 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)\text d\varphi\text d\theta \end{align*}

Um nun die Integrale ausrechnen zu können, wird die explizite Form des Spannungstensors für dieses Problem benötigt. Im Allgemeinen ist der Spannungstensor für reale, newtonsche Fluide gegeben durch: \begin{align*} \underline{\underline\sigma}= -\bigg[p-\zeta(\vec{\nabla}\cdot\vec{u})\bigg]\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*}

Mit der Volumenviskosität \(\zeta\), welche jedoch aufgrund der Stokes'schen Hypothese vernachlässigt werden kann. Damit folgt: \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*}

Dabei ist \(\mathbb{1}\) der Einheitstensor, \(\mu\) die dynamische Viskosität und \(\vec{u}\) die Geschwindigkeit des Fluides. Aufgrund dessen, dass bei der Kraftberechnung in Kugelkoordinaten nur die \(rr\)-, \(\varphi r\)- und \(\theta r\)-Komponente des Spannungstensors eingehen, werden im Folgenden auch nur diese bestimmt. 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(R,\theta) +\mu(R,\theta)\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(R,\theta) +\mu(R,\theta)\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*}

Die Verwendung der No-Slip-RB bedeutet, dass sowohl die normale als auch die tangentiale Komponente des Geschwindigkeitsvektors auf der Oberfläche eines festen Körpers verschwinden müssen. Damit folgt: \begin{align*} u_r\Big|_{r=R}=u_{\theta}\Big|_{r=R}=0 \end{align*}

Außerdem impliziert das Verschwinden des Geschwindigkeitsvektors 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(R,\theta) +\frac{4}{3}\mu(R,\theta) \frac{\partial u_r}{\partial r}\Big|_{r=R} \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*} \sigma_{\theta r}\Big|_{r=R}= \mu(R,\theta)\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*}

Mit den Randbedingungen folgt dann final für die \(\theta r\)-Komponente: \begin{align*} \sigma_{\theta r}\Big|_{r=R}= \mu(R,\theta) \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(R,\theta) +\frac{4}{3}\mu(R,\theta)\frac{\partial u_r}{\partial r}\Big|_{r=R}\bigg)\vec{e}_r +\mu(R,\theta) \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(R,\theta)\vec{e}_r\sin(\theta)d\theta +\frac{4}{3}\int_{0}^{\pi}\mu(R,\theta) \frac{\partial u_r}{\partial r}\Big|_{r=R} \vec{e}_r\sin(\theta)\text d\theta +\int_{0}^{\pi}\mu(R,\theta) \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(R,\theta)\vec{e}_r\cdot\vec{e}_z \sin(\theta)d\theta +\frac{4}{3}\int_{0}^{\pi}\mu(R,\theta) \frac{\partial u_r}{\partial r}\Big|_{r=R} \vec{e}_r\cdot\vec{e}_z\sin(\theta)\text d\theta +\int_{0}^{\pi}\mu(R,\theta) \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(R,\theta)\cos(\theta)\sin(\theta)d\theta +\frac{4}{3}\int_{0}^{\pi}\mu(R,\theta) \frac{\partial u_r}{\partial r}\Big|_{r=R} \cos(\theta)\sin(\theta)\text d\theta -\int_{0}^{\pi}\mu(R,\theta) \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 +\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*}

Natürliche Einheiten

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.

Viskosität, dynamisch vs. kinematisch

Von den beiden, \(\mu=\nu\rho\), ist keine ein wirklicher Materialparameter, doch die dynamische Viskosität hängt nur von der Temperatur (und von WW-Details der Fluid-Partikel) ab. Da im Stokes-Regime die Temperatur-Änderung klein sein sollte, ist es also sinnvoller, \(\mu\) als "Materialparameter" vorzugeben.

Die Viskosität wird in der Funktion Visc_nu() in visc_nu.c berechnet, ein Aufruf pro Zelle. Trotz der Namen nu1 (Scherv.) und nu2 (Volumenv.) handelt es sich um die dynamischen Viskositäten. Beim Vorgeben von \(\mu\) sollte \(\rho\) (im Code v[RHO]) also gar nicht auftauchen. Verwendet man obige "zähe Einheiten", so wird die Funktion besonders einfach:

void Visc_nu(double *v, double x1, double x2, double x3, double *nu1, double *nu2)
/* ... */
{
  *nu1 = 1.0;  *nu2 = 0.0;
}

Analytische Lösung

Mit \(R=1\) und \(\vec u_\infty=-v\vec e_z\):

$$ \frac{\vec u}{v} = \frac{\cos\theta}{2}(-r^{-3}+3r^{-1}-2)\vec e_r+\frac{\sin\theta}{4}(-r^{-3}-3r^{-1}+4)\vec e_\theta $$

Entgegen "der" Intuition resultiert daraus $$ \vec\nabla\times\vec u=\frac{3v}{2r^2}\sin\theta\,\vec e_\varphi~. $$

Jenseits von Stokes: Oseen-Fluss

Implementierung in Python

Nachdem die Daten mittels der Bibliothek pyPLUTO geladen wurden, wurde zunächst ein Meshgrid für die Radien und die Winkel erstellt:

RS, THETAS = np.meshgrid(D.x1, D.x2)

Hierbei wurde verwendet:

D = pp.pload(500, w_dir="C:/Users/...")

Die oben beschriebenen Randbedingungen sind jedoch nicht in den Daten für Radien und den Winkeln bzw. für die jeweiligen Geschwindigkeitskomponenten enthalten. Daher müssen sie manuell in die Datensätze eingefügt werden. Dafür sorgt folgende Funktion:

def before_first_value_line(array, value):
    
    new_arr = np.hstack((np.full((array.shape[0], 1), value), array))
    
    return new_arr

np.hstack() ist dabei eine Funktion, um zwei Arrays nebeneinander zu setzen, sodass zwei Arrays der Form \((n,m)\) und \((n,p)\) zu einem Array der Form \((n,m+p)\) zusammengefügt werden. np.full() erstellt hierbei ein Array mit der selben Anzahl an Zeilen, wie das zu übergebende Array, jedoch nur mit einer Spalte.

Mit dieser Funktion lassen sich jetzt die Randbedingungen in die Datensätze einfügen:

R = before_first_value_line(RS,1) 
THETA = THETAS

VR = before_first_value_line(D.vx1.T,0) 
VTHETA = before_first_value_line(D.vx2.T,0)

Um nun \(F_z\) gemäß der obigen Gleichung bestimmen zu können, müssen zunächst die partiellen Ableitungen der Geschwindigkeitskomponenten nach Radius und Polarwinkel numerisch bestimmt werden. Dafür erweist sich die Numpy-Funktion np.gradient() als nützlich. Mit dieser lässt sich eine Funktion definieren, die die Ableitung nach dem Radius durchführt:

def partial_derivative_dr(V, R, i, j):
    
    dvdr = np.gradient(V[i,:], R[i,:])
    
    return dvdr[j]

Dabei durchläuft der Index \(i\) die Polarwinkel und der Index \(j\) die Radien. Zusätzlich muss jedoch darauf geachtet werden, dass die Transponierten der Geschwindigkeits-Arrays für \(V\) eingesetzt werden müssen. Um nun die Kraft \(dF_z\) in \(z\)-Richtung auf eine Zelle gemäß obiger Gleichung bestimmen zu können müsste zusätzlich die dynamische Viskosität \(\mu\) bestimmt werden. Werden jedoch zähe Einheiten verwendet, so ist \(\mu=1\). Damit lässt sich nun eine Funktion schreiben, die die Kraft \(dF_z\) in \(z\)-Richtung auf eine Zelle gemäß obiger Gleichung bestimmt:

def local_force(VR, VTHETA, i, j, Ma, Re, R, P):
    
    df = - P[i,j] * np.cos(THETA[i,j]) * np.sin(THETA[i,j]) * np.pi / (len(THETA[:,0]))
    df += dynamic_viscosity(Ma,Re,R) * (4/3) * partial_derivative_dr(VR,R,i,j) * np.cos(THETA[i,j]) * np.sin(THETA[i,j]) * np.pi / (len(THETA[:,0]))
    df += dynamic_viscosity(Ma,Re,R) * partial_derivative_dr(VTHETA,R,i,j) * np.sin(THETA[i,j]) ** 2 * np.pi / (len(THETA[:,0]))
    
    return df

Dabei ist \(P\) der Durck, welcher von der jeweiligen Zelle aus auf die Kugel ausgeübt wird. Dieser lässt sich ebenfalls aus den Daten laden. Der Faktor np.pi / (len(THETA[:,0])) entspricht hierbei dem \(d\theta\) in obiger Gleichung. Für die Gesamtkraft \(F_z\) in \(z\)-Richtung auf die Kugel müssen die \(dF_z\) für alle Zellen um die Kugel aufsummiert werden und anschließend mit dem Faktor \(2\pi R^2\) gewichtet werden. Außerdem muss das Endergebnis zusätzlich mit einem Faktor \(2\) multipliziert werden, da aus Symmetriegründen nur die Hälfte der Kugel simuliert wird. Damit lässt sich folgende Funktion definieren:

def net_force(VR, VTHETA, j, Ma, Re, R, P): 
    
    df_arr = np.zeros(len(THETA[:,0])) 
    
    for i in np.arange(0,len(THETA[:,0])):
        
        df_arr[i] = local_force(VR, VTHETA, i, j, Ma, Re, R, P)
        
    f = 2 * np.pi * np.min(R) ** 2 * 2 * np.sum(df_arr) 

    return f