[PYTHON] Numerische Berechnung von kompressiblem Fluid nach der Methode des endlichen Volumens

1. Zuallererst

Die Kompressibilität von Flüssigkeiten muss für Strömungen um sich schnell bewegende Objekte wie Flugzeuge und Raketen sowie für Phänomene mit starken Druck-, Temperatur- und Dichteänderungen wie Verbrennung und Explosion berücksichtigt werden. Kompressible Flüssigkeiten verursachen verschiedene Phänomene wie Schallwellen, Stoßwellen und aerodynamische Erwärmung, die in inkompressiblen Flüssigkeiten nicht zu finden sind, und ihre numerischen Berechnungen haben andere Schwierigkeiten als inkompressible Flüssigkeiten. In diesem Artikel möchte ich eine Methode zur numerischen Berechnung von komprimierbarem Fluid mit der Finite-Volumen-Methode (FVM) vorstellen.

2. Formulierung der Methode des endlichen Volumens

Betrachten Sie das Folgende als partielle Differentialgleichung für den Vektor $ \ mathbf {U} $, der ein Array von $ N $ unabhängigen Variablen ist.

\frac{\partial \mathbf{U}}{\partial t} +
\frac{\partial \mathbf{F}}{\partial x} +
\frac{\partial \mathbf{G}}{\partial y} +
\frac{\partial \mathbf{H}}{\partial z} = \mathbf{S}

Zu diesem Zeitpunkt heißt $ \ mathbf {F}, \ mathbf {G}, \ mathbf {H} $ ** Flux **, und $ \ mathbf {S} $ repräsentiert den Quellbegriff von außen. Erwägen Sie, diese Gleichung numerisch nach der Methode des endlichen Volumens zu lösen. Die Methode des endlichen Volumens unterteilt den Berechnungsbereich in Zellen, die als ** Kontrollvolumina ** bezeichnet werden. Dieses Kontrollvolumen muss kein gleichmäßiges Quadrat sein, und unstrukturierte Gitter können verwendet werden. Wie bei der Finite-Elemente-Methode besteht eine der großen Stärken der Finite-Volumen-Methode darin, dass sie komplexe Formen verarbeiten kann, die realen Problemen entsprechen. Nun sei die obige Gleichung das $ n $ -te Mal $ t_n $ zum $ n + 1 $ -ten Mal $ t_ {n + 1} $ (vorausgesetzt, der Zeitschritt ist $ \ Delta t $) in Bezug auf die Zeit. Integrieren Sie in das $ i $ -te Kontrollvolumen $ V_i $ für den Speicherplatz (Referenz [^ 1], [^ 2]).

\int_{t_n}^{t_{n+1}}dt\iiint_{V_i}dV\frac{\partial \mathbf{U}}{\partial t} + \int_{t_n}^{t_{n+1}}dt\iiint_{V_i}dV\left[\frac{\partial \mathbf{F}}{\partial x} + 
\frac{\partial \mathbf{G}}{\partial y} + 
\frac{\partial \mathbf{H}}{\partial z}\right] = \int_{t_n}^{t_{n+1}}dt\iiint_{V_i}dV\mathbf{S}

Die zeitliche Integration eines Elements auf der linken Seite kann sofort durchgeführt werden, und unter Verwendung des Gaußschen Divergenzsatzes wird die Volumenintegration von zwei Elementen auf der linken Seite in den Bereich auf der Oberfläche $ \ partiell V_i $ des Kontrollvolumens umgewandelt.

\iiint_{V_i}dV\left[\mathbf{U}^{n+1}-\mathbf{U}^n\right] + \int_{t_n}^{t_{n+1}}dt\iint_{\partial V_i}d\Omega\left[\mathbf{F}n_x + \mathbf{G}n_y + \mathbf{H}n_z\right] =
\int_{t_n}^{t_{n+1}}dt\iiint_{V_i}dV\mathbf{S}

$ \ Mathbf {U} ^ n = \ mathbf {U} \ left (t_n, \ mathbf {x} \ right) $. Wenn diese Fläche $ \ partiell V_i $ aus mehreren weiteren Flächen $ \ Omega_ {i, l} $ besteht (z. B. 6 Rechtecke für ein Quadrat), schreiben Sie die Fläche in die Summe jeder Fläche. Ich werde.

\iiint_{V_i}dV\left[\mathbf{U}^{n+1}-\mathbf{U}^n\right] + \sum_{l}\int_{t_n}^{t_{n+1}}dt\iint_{\Omega_{i,l}}d\Omega\left[\mathbf{F}n_x + \mathbf{G}n_y + \mathbf{H}n_z\right] =
\int_{t_n}^{t_{n+1}}dt\iiint_{V_i}dV\mathbf{S}

Deshalb,

\mathbf{U}_i^n = \frac{1}{\Delta V_i}\iiint_{V_i}dV\mathbf{U}^n,\ 
\mathbf{S}_i^n = \frac{1}{\Delta t\Delta V_i}\int_{t_n}^{t_{n+1}}dt\iiint_{V_i}dV\mathbf{S}

\mathbf{F}_{i,l}^n = \frac{1}{\Delta t\Delta \Omega_{i,l}}\int_{t_n}^{t_{n+1}}dt\iint_{\Omega_{i,l}}d\Omega\mathbf{F},\ 
\mathbf{G}_{i,l}^n = \frac{1}{\Delta t\Delta \Omega_{i,l}}\int_{t_n}^{t_{n+1}}dt\iint_{\Omega_{i,l}}d\Omega\mathbf{G},\ 
\mathbf{H}_{i,l}^n = \frac{1}{\Delta t\Delta \Omega_{i,l}}\int_{t_n}^{t_{n+1}}dt\iint_{\Omega_{i,l}}d\Omega\mathbf{H}

Wenn definiert als, erhalten wir:

\frac{\mathbf{U}_i^{n+1}-\mathbf{U}_i^n}{\Delta t} + \sum_{l}\left[\mathbf{F}_{i,l}^n n_x+\mathbf{G}_{i,l}^n n_y+\mathbf{H}_{i,l}^n n_z\right]\frac{\Delta \Omega_{i,l}}{\Delta V_i} = \mathbf{S}_i^n

Zu diesem Zeitpunkt ist $ \ mathbf {F} ^ n_ {i, l}, \ \ mathbf {G} ^ n_ {i, l}, \ \ mathbf {H} ^ n_ {i, l} $ ** ein numerischer Fluss Es wird ein Bündel ** genannt. Dies ist die häufigste Form der Methode mit endlichem Volumen. Wenn das Kontrollvolumen gleichmäßig in Quadrate $ V_ {i, j, k} $ auf jeder Seite von $ \ Delta x, \ Delta y, \ Delta z $ unterteilt ist, kann es wie folgt transformiert werden.

\frac{\mathbf{U}_{i,j,k}^{n+1}-\mathbf{U}_{i,j,k}^n}{\Delta t} +
\frac{\mathbf{F}_{i+1/2,j,k}^n-\mathbf{F}_{i-1/2,j,k}^n}{\Delta x} +
\frac{\mathbf{G}_{i,j+1/2,k}^n-\mathbf{G}_{i,j-1/2,k}^n}{\Delta y} + 
\frac{\mathbf{H}_{i,j,k+1/2}^n-\mathbf{H}_{i,j,k-1/2}^n}{\Delta z} = \mathbf{S}_{i,j,k}^n

Es ist ziemlich vertraut geworden. Die Indizes $ \ mathbf {F} ^ n_ {i + 1/2, j, k} $ hier sind $ V_ {i, j, k} $ und $ V_ {i + 1, j, k} $ Dies bedeutet, dass der Wert der Grenze verwendet wird. Alle oben genannten sind mathematisch exakte Varianten, aber in der Realität ist es nicht möglich, jeden Fluss auf der linken Seite genau zu bestimmen. Das Wesentliche der Methode des endlichen Volumens ist daher die Frage, wie der numerische Fluss an dieser Grenzfläche angenähert werden kann. Unten sehen Sie eine schematische Darstellung des Falls, in dem jedes Reglervolumen in zwei Dimensionen rechteckig ist. Der Wert im Kontrollvolumen wird durch das "Gleichgewicht" der Menge bestimmt, die durch die Grenzfläche ein- und ausgeht.

3. Eigenschaften der komprimierbaren Flüssigkeit

3-1. Gleichungssystem für kompressible Flüssigkeit

Lassen Sie uns hier das Gleichungssystem für komprimierbare Flüssigkeiten einmal organisieren. Unter Berücksichtigung jedes Erhaltungsgesetzes von Dichte, Impuls und Energie kann die Gleichung in der am Anfang von 2 angegebenen Form wie folgt geschrieben werden (eine detaillierte Ableitung finden Sie im Lehrbuch der Fluiddynamik).

\mathbf{U} = \begin{bmatrix}
\rho \\
\rho u \\
\rho v \\
\rho w \\
\rho E
\end{bmatrix},\ 
\mathbf{F} = \begin{bmatrix}
\rho u \\
\rho u^2+p \\
\rho uv \\
\rho uw \\
\rho uH
\end{bmatrix},\ 
\mathbf{G} = \begin{bmatrix}
\rho v \\
\rho vu \\
\rho v^2+p \\
\rho vw \\
\rho vH
\end{bmatrix},\ 
\mathbf{H} = \begin{bmatrix}
\rho w \\
\rho wu \\
\rho wv \\
\rho w^2+p \\
\rho wH
\end{bmatrix}

E = \frac{1}{\gamma-1}\frac{p}{\rho}+\frac{1}{2}\left(u^2+v^2+w^2\right),\ H = E+\frac{p}{\rho}

Hier sind $ \ rho und p $ Dichte und Druck, $ u, v und w $ sind Geschwindigkeiten in $ x-, y- und z $ -Richtung und $ E und H $ sind Energie und Enthalpie pro Masseneinheit. Repräsentiert. $ \ Gamma $ ist das spezifische Wärmeverhältnis. Auch die Begriffe Viskosität / Wärmeleitung und Quelle werden hier ignoriert.

3-2. Jacobi-Diagonalisierung

Ein wichtiger Unterschied zwischen komprimierbaren und nicht komprimierbaren Flüssigkeiten ist das Vorhandensein von Druckwellen oder Schallwellen. Daher werden die Informationen zusätzlich zu den entlang des Flusses übertragenen Informationen auch durch Schallwellen übertragen. Für jeden Fluss ist also der Jacobi $ \ mathbf {A} = \ partiell \ mathbf {F} / \ partiell \ mathbf {U}, \ mathbf {B} = \ partiell \ mathbf {G} / \ partiell \ mathbf Versuchen Sie, {U}, \ mathbf {C} = \ teilweise \ mathbf {H} / \ teilweise \ mathbf {U} $ zu finden. Für $ \ mathbf {A} $ wird es beispielsweise wie folgt berechnet (Referenz [^ 3]).

\mathbf{A} = \begin{bmatrix}
0 & 1 & 0 & 0 & 0 \\
\left(\gamma-1\right)H-u^2-a^2 & -\left(\gamma-3\right)u & -\left(\gamma-1\right)v & -\left(\gamma-1\right)w & \gamma-1 \\
-uv & v & u & 0 & 0 \\
-uw & w & 0 & u & 0 \\
\frac{1}{2}u\left[\left(\gamma-3\right)H-a^2\right] & H-\left(\gamma-1\right)u^2 & -\left(\gamma-1\right)uv & -\left(\gamma-1\right)uw & \gamma u
\end{bmatrix}

Wobei $ a = \ sqrt {\ gamma p / \ rho} $, was die Schallgeschwindigkeit darstellt. Alle diese Jacobi-Eigenwerte sind real und können mit der folgenden Matrix $ \ mathbf {K} $ diagonalisiert werden (diese Berechnung ist sehr mühsam, ich möchte sie nicht noch einmal durchführen).

\mathbf{K} = \begin{bmatrix}
   1 &                                   1 & 0 & 0 &    1 \\
 u-a &                                   u & 0 & 0 &  u+a \\
   v &                                   v & 1 & 0 &    v \\
   w &                                   w & 0 & 1 &    w \\
H-ua & \frac{1}{2}\left(u^2+v^2+w^2\right) & v & w & H+ua
\end{bmatrix}\\
\mathbf{K}^{-1} = \frac{\gamma-1}{2a^2}\begin{bmatrix}
H+\frac{a}{\gamma-1}\left(u-a\right) & -u-\frac{a}{\gamma-1} & -v & -w & 1 \\
-2H+\frac{4a^2}{\gamma-1} & 2u & 2v & 2w & -2 \\
-\frac{2va^2}{\gamma-1} & 0 & \frac{2a^2}{\gamma-1} & 0 & 0 \\
-\frac{2wa^2}{\gamma-1} & 0 & 0 & \frac{2a^2}{\gamma-1} & 0 \\
H-\frac{a}{\gamma-1}\left(u+a\right) & -u+\frac{a}{\gamma-1} & -v & -w & 1
\end{bmatrix}\\
\mathbf{\Lambda} = \begin{bmatrix}
u-a & 0 & 0 & 0 &   0 \\
  0 & u & 0 & 0 &   0 \\
  0 & 0 & u & 0 &   0 \\
  0 & 0 & 0 & u &   0 \\
  0 & 0 & 0 & 0 & u+a
\end{bmatrix}\\
\mathbf{\Lambda} = \mathbf{K}^{-1}\mathbf{A}\mathbf{K}, \mathbf{A} = \mathbf{K}\mathbf{\Lambda}\mathbf{K}^{-1}

3-3. Charakteristische Geschwindigkeit und Lehman-Invariante

Betrachten wir der Einfachheit halber den eindimensionalen Fall. Das Gleichungssystem für eine Dimension ist wie folgt.

\frac{\partial \mathbf{U}}{\partial t} +
\frac{\partial \mathbf{F}}{\partial x} = \mathbf{0}

\mathbf{U} = \begin{bmatrix}
\rho \\
\rho u \\
\rho E
\end{bmatrix},\ 
\mathbf{F} = \begin{bmatrix}
\rho u \\
\rho u^2+p \\
\rho uH
\end{bmatrix}

E = \frac{1}{\gamma-1}\frac{p}{\rho}+\frac{1}{2}u^2,\ H = E+\frac{p}{\rho}

Zu diesem Zeitpunkt können wir durch Einführen der folgenden Beträge $ I_0, I_ \ pm $ zeigen, dass die folgende Gleichung gilt.

I_0 = \frac{p}{\rho^\gamma},\ 
I_\pm = u \pm \int\frac{a}{\rho}d\rho

\frac{\partial I_0}{\partial t} + u\frac{\partial I_0}{\partial x} = 0 \\
\frac{\partial I_\pm}{\partial t} + \left(u \pm a\right)\frac{\partial I_\pm}{\partial x} = 0

Unter der Annahme, dass $ I_0 $ immer konstant ist, dh die Bedingung ** gleiche Entropie **, wird $ I_ \ pm $ wie folgt vereinfacht.

I_\pm = u \pm \frac{2a}{\gamma-1}

Dies wird als ** Lehman-Invariante ** bezeichnet, was bedeutet, dass sich $ I_0, I_ + und I_- $ mit Geschwindigkeiten $ u, u + a bzw. u-a $ ausbreiten. Mit anderen Worten wurde gefunden, dass die wahre Natur der in 3-2 erhaltenen Eigenwerte die Ausbreitungsgeschwindigkeit dieser Invarianten war. Daher wird dieser Eigenwert als ** charakteristische Geschwindigkeit ** bezeichnet, da es die Geschwindigkeit ist, die das Gleichungssystem charakterisiert. Bei einer Dimension gibt es drei unabhängige Variablen. Wenn also $ I_0, I_ +, I_- $ gefunden werden, können alle Variablen gefunden werden. Mit anderen Worten kann gesagt werden, dass der Wert jedes Punktes durch die Kombination der drei Arten von Wellen bestimmt wird, die sich ausgebreitet haben. Die folgende Abbildung zeigt dies schematisch.

4. Anwendung auf die Methode des endlichen Volumens

4-1. Gegenwindunterschied

Betrachten Sie zur Vereinfachung die folgende eindimensionale skalare Übertragungsgleichung.

\frac{\partial u}{\partial t} + \frac{\partial f}{\partial x} = 0,\ f = cu

Wenn Sie dies nach der Methode des endlichen Volumens dispergieren, erhalten Sie Folgendes.

\frac{u_i^{n+1}-u_i^n}{\Delta t} + \frac{f_{i+1/2}^n-f_{i-1/2}^n}{\Delta x} = 0

Zu diesem Zeitpunkt repräsentiert $ c $ die Übertragungsgeschwindigkeit, und die Methode zum Zuweisen des Flusswerts gemäß dem Vorzeichen wird als ** Aufwinddifferenzmethode ** bezeichnet.

f_{i+1/2}^n = \begin{cases}
cu_{i}^n & (c \ge 0) \\
cu_{i+1}^n & (c < 0)
\end{cases}

In der Translokationsgleichung wird der Wert immer durch die Informationen aus dem Upstream bestimmt, sodass der Wert auf diese Weise ausgewählt wird. Wie können wir dies nun für mehrere Variablen erweitern? Wie bereits erwähnt, gibt es bei mehreren Variablen mehrere Arten von Wellen, daher muss diese gut getrennt werden. Da es möglich ist, dass sich einige Wellen vorwärts und einige Wellen rückwärts ausbreiten, ist es notwendig, die Richtung der Aufwinddifferenz in Abhängigkeit von der Welle zu ändern. Daher werden wir einen neuen Vektor $ \ mathbf {W} = \ mathbf {K} ^ {-1} \ mathbf {U} $ unter Verwendung der in 3-2 berücksichtigten Diagonalisierung einführen. Unter der Annahme, dass die Änderung in $ \ mathbf {K} $ gering ist, kann die Gleichung wie folgt transformiert werden.

\frac{\partial\mathbf{W}}{\partial t} + \mathbf{\Lambda}\frac{\partial\mathbf{W}}{\partial x} = \mathbf{0}

Daher sind 5 echte eindeutige Werte $ \ lambda_1 = ua, \ lambda_2 = \ lambda_3 = \ lambda_4 = u, \ lambda_5 = u + a $ und $ \ mathbf {W} = \ left [W_1, W_2, W_3, W_4, W_5 \ right] $ transformiert die Gleichung in fünf unabhängige Translokationsgleichungen:

\frac{\partial W_l}{\partial t} + \lambda_l \frac{\partial W_l}{\partial x} = 0

Dann scheint es gut, die Aufwinddifferenz separat auf jede dieser fünf Translokationsgleichungen anzuwenden. Wenn Sie also den Betrag $ \ lambda_l ^ + = \ max \ left (\ lambda_l, 0 \ right), \ lambda_l ^ - = \ min \ left (\ lambda_l, 0 \ right) $ einführen, wird diese Translokationsgleichung gewickelt. Die durch die obige Differenz diskreten Gleichungen können wie folgt zu einer Gleichung kombiniert werden.

\frac{\left(W_l\right)_i^{n+1}-\left(W_l\right)_i^n}{\Delta t}
+ \lambda_l^+\frac{\left(W_l\right)_i^n-\left(W_l\right)_{i-1}^n}{\Delta x}
+ \lambda_l^-\frac{\left(W_l\right)_{i+1}^n-\left(W_l\right)_i^n}{\Delta x} = 0

Dann $ \ mathbf {\ Lambda} ^ + = \ mathrm {diag} \ left (\ lambda_1 ^ +, \ lambda_2 ^ +, \ lambda_3 ^ +, \ lambda_4 ^ +, \ lambda_5 ^ + \ right), \ mathbf {\ Lambda} ^ - = \ mathrm {diag} \ left (\ lambda_1 ^ -, \ lambda_2 ^ -, \ lambda_3 ^ -, \ lambda_4 ^ -, \ lambda_5 ^ - \ right) $ Schreiben Sie in Form.

\frac{\mathbf{W}_i^{n+1}-\mathbf{W}_i^n}{\Delta t}
+ \mathbf{\Lambda}^+ \frac{\mathbf{W}_i^n-\mathbf{W}_{i-1}^n}{\Delta x}
+ \mathbf{\Lambda}^- \frac{\mathbf{W}_{i+1}^n-\mathbf{W}_i^n}{\Delta x} = \mathbf{0}

Schließlich $ \ mathbf {A} ^ + = \ mathbf {K} \ mathbf {\ Lambda} ^ + \ mathbf {K} ^ {-1}, \ mathbf {A} ^ - = \ mathbf {K} \ Als mathbf {\ Lambda} ^ - \ mathbf {K} ^ {-1} $ lautet die Formel für $ \ mathbf {U} $ wie folgt (unter der Annahme, dass die Änderung in Jacobian gering ist). Der Index der Position wird getäuscht)

\frac{\mathbf{U}_i^{n+1}-\mathbf{U}_i^n}{\Delta t}
+ \mathbf{A}^+ \frac{\mathbf{U}_i^n-\mathbf{U}_{i-1}^n}{\Delta x}
+ \mathbf{A}^- \frac{\mathbf{U}_{i+1}^n-\mathbf{U}_i^n}{\Delta x} = \mathbf{0}

4-2. Methode zur Trennung der Flussflussdifferenz

$\mathbf{A}^+ \left(\mathbf{U}^n_i-\mathbf{U}^n_{i-1}\right) = \Delta{\mathbf{F}^+}^n_i,\ \mathbf{A}^- \left(\mathbf{U}^n_{i+1}-\mathbf{U}^n_i\right) = \Delta{\mathbf{F}^-}^n_i,\ Wenn Sie \ Delta {\ mathbf {F} ^ +} ^ n_i + \ Delta {\ mathbf {F} ^ -} ^ n_i = \ Delta \ mathbf {F} ^ n_i $ schreiben, lautet der letzte Ausdruck in 4-1

\frac{\mathbf{U}_i^{n+1}-\mathbf{U}_i^n}{\Delta t} + \frac{\Delta\mathbf{F}_i^n}{\Delta x} =
\frac{\mathbf{U}_i^{n+1}-\mathbf{U}_i^n}{\Delta t} + \frac{\Delta{\mathbf{F}^+}_i^n}{\Delta x} + \frac{\Delta{\mathbf{F}^-}_i^n}{\Delta x} = \mathbf{0}

Es ist ersichtlich, dass die Flussflussdifferenz $ \ Delta \ mathbf {F} ^ n_i $ in positiver und negativer Richtung getrennt ist. Die auf dieser Formel basierende Berechnungsmethode heißt ** Flux Difference Splitting (FDS) **. Lassen Sie uns dies ein wenig mehr transformieren.\mathbf{A}^++\mathbf{A}^- = \mathbf{A}Und auch\mathbf{A}^+-\mathbf{A}^- = \left|\mathbf{A}\right|Wenn definiert als(Der Index n wird für eine Weile weggelassen)、

\mathbf{A}^+ \frac{\mathbf{U}_i-\mathbf{U}_{i-1}}{\Delta x} + \mathbf{A}^- \frac{\mathbf{U}_{i+1}-\mathbf{U}_i}{\Delta x}
= \frac{1/2\mathbf{A}\left(\mathbf{U}_{i+1} + \mathbf{U}_i\right)-1/2\left|\mathbf{A}\right|\left(\mathbf{U}_{i+1}-\mathbf{U}_i\right)-1/2\mathbf{A}\left(\mathbf{U}_i+\mathbf{U}_{i-1}\right)+1/2\left|\mathbf{A}\right|\left(\mathbf{U}_i-\mathbf{U}_{i-1}\right)}{\Delta x}

Also der Fluss

\begin{eqnarray}
\mathbf{F}_{i+1/2}
&=& \frac{1}{2}\mathbf{A}\left(\mathbf{U}_{i+1} + \mathbf{U}_i\right)-\frac{1}{2}\left|\mathbf{A}\right|\left(\mathbf{U}_{i+1}-\mathbf{U}_i\right)\\
&=& \frac{1}{2}\left(\mathbf{F}\left(\mathbf{U}_{i+1}\right) + \mathbf{F}\left(\mathbf{U}_i\right)\right)-\frac{1}{2}\left|\mathbf{A}\right|\left(\mathbf{U}_{i+1}-\mathbf{U}_i\right)
\end{eqnarray}

Kann geschrieben werden als. Übrigens habe ich die ganze Zeit den Index von Jacobians Position getäuscht, aber wie soll ich ihn tatsächlich berechnen? Daher führte Roe einen speziellen Durchschnitt ein, der wie folgt mit der Quadratwurzel der Dichte gewichtet ist [^ 4].

\begin{eqnarray}
\bar{\rho}_{i+1/2} &=& \sqrt{\rho_{i+1}\rho_i} \\
\bar{u}_{i+1/2} &=& \frac{\sqrt{\rho_{i+1}}u_{i+1}+\sqrt{\rho_i}u_i}{\sqrt{\rho_{i+1}}+\sqrt{\rho_i}}\\
\bar{v}_{i+1/2} &=& \frac{\sqrt{\rho_{i+1}}v_{i+1}+\sqrt{\rho_i}v_i}{\sqrt{\rho_{i+1}}+\sqrt{\rho_i}}\\
\bar{w}_{i+1/2} &=& \frac{\sqrt{\rho_{i+1}}w_{i+1}+\sqrt{\rho_i}w_i}{\sqrt{\rho_{i+1}}+\sqrt{\rho_i}}\\
\bar{H}_{i+1/2} &=& \frac{\sqrt{\rho_{i+1}}H_{i+1}+\sqrt{\rho_i}H_i}{\sqrt{\rho_{i+1}}+\sqrt{\rho_i}}\\
\bar{a}_{i+1/2} &=& \sqrt{\left(\gamma-1\right)\left[\bar{H}_{i+1/2}-\frac{1}{2}\left(\bar{u}_{i+1/2}^2+\bar{v}_{i+1/2}^2+\bar{w}_{i+1/2}^2\right)\right]}
\end{eqnarray}

Mit ihnen\left|\bar{\mathbf{A}}_{i+1/2}\right|Und ersetzen Sie den folgenden Fluss.

\mathbf{F}_{i+1/2}
= \frac{1}{2}\left(\mathbf{F}\left(\mathbf{U}_{i+1}\right) + \mathbf{F}\left(\mathbf{U}_i\right)\right)-\frac{1}{2}\left|\bar{\mathbf{A}}_{i+1/2}\right|\left(\mathbf{U}_{i+1}-\mathbf{U}_i\right)

Dies wird als ** Roes ungefähre Lehman-Lösung ** bezeichnet.

4-3. Strömungsfluss-Trennmethode

Als andere Methode wird der Fluss direkt getrennt, ohne den Jacobi wie folgt zu durchlaufen, und den Jacobi $ \ partiell \ mathbf {F} ^ + / \ partiell \ mathbf {U}, \ partiell \ mathbf {F} ^ Es ist denkbar, die eindeutigen Werte von- / \ partiell \ mathbf {U} $ so zu konfigurieren, dass sie nicht negativ bzw. nicht positiv sind.

\frac{\mathbf{U}_i^{n+1}-\mathbf{U}_i^n}{\Delta t} +
\frac{\mathbf{F}^+\left(\mathbf{U}_i^n\right)-\mathbf{F}^+\left(\mathbf{U}_{i-1}^n\right)}{\Delta x} +
\frac{\mathbf{F}^-\left(\mathbf{U}_{i+1}^n\right)-\mathbf{F}^-\left(\mathbf{U}_i^n\right)}{\Delta x} = \mathbf{0}

\mathbf{F}^+\left(\mathbf{U}_i^n\right)+\mathbf{F}^-\left(\mathbf{U}_{i+1}^n\right) = \mathbf{F}_{i+1/2}^n

Diese Art von Methode wird als ** Flux Vector Splitting (FVS) ** bezeichnet.

4-3-1. Steger-Erwärmungsschema

Die einfachste Methode besteht darin, den Fluss wie folgt zu trennen, wobei auf die Formel in 4-1 Bezug genommen wird.

\mathbf{F}_{i+1/2}^n = \mathbf{F}^+\left(\mathbf{U}_i^n\right)+\mathbf{F}^-\left(\mathbf{U}_{i+1}^n\right) = \mathbf{A}^+\left(\mathbf{U}_i^n\right)\mathbf{U}_i^n+\mathbf{A}^-\left(\mathbf{U}_{i+1}^n\right)\mathbf{U}_{i+1}^n

Dies wird als ** Steger-Warming-Schema ** [^ 5] bezeichnet. Darüber hinaus wurden verschiedene Verfahren wie das Strömungsfluss-Trennverfahren von van Leer [^ 6] vorgeschlagen.

4-3-2. AUSM Fluss in Richtung $ x $

\mathbf{F} = \begin{bmatrix}
\rho u \\
\rho u^2+p \\
\rho uv \\
\rho uw \\
\rho uH
\end{bmatrix}

Wenn Sie sich das ansehen, können Sie sehen, dass alles außer Druck eine Übertragungsform hat, die proportional zur Strömungsgeschwindigkeit $ u $ ist. Daher schlugen Liou und Steffen eine Methode namens ** Advance Upstream Splitting Method (AUSM) ** vor, die den Fluss wie folgt in eine Übertragungskomponente und eine Druckkomponente unterteilt [^ 7].

\mathbf{F}^{\left(a\right)} = u\begin{bmatrix}
\rho \\
\rho u \\
\rho v \\
\rho w \\
\rho H
\end{bmatrix} = M\begin{bmatrix}
\rho a \\
\rho au \\
\rho av \\
\rho aw \\
\rho aH
\end{bmatrix},\ 
\mathbf{F}^{\left(p\right)} = \begin{bmatrix}
0 \\
p \\
0 \\
0 \\
0
\end{bmatrix}

\mathbf{F}^{\left(a\right)}+\mathbf{F}^{\left(p\right)} = \mathbf{F}

Hier heißt $ M = u / a $ ** Machzahl ** und ist ein wichtiger Parameter bei der Betrachtung komprimierbarer Flüssigkeiten. In AUSM wird der Absolutwert dieser Machzahl $ M $ so eingestellt, dass der Ausdruck an der Grenze von 1 umgeschaltet wird (dh ob es sich um Überschall handelt oder nicht). Führen Sie zunächst für die Übertragungskomponente $ \ mathbf {F} ^ {\ left (a \ right)} $ die folgenden Schritte aus (der Index n wird für eine Weile weggelassen).

{\mathbf{F}^{\left(a\right)}}_{i+1/2} = \max\left(M_{i+1/2}, 0\right)\begin{bmatrix}
\rho a \\
\rho au \\
\rho av \\
\rho aw \\
\rho aH
\end{bmatrix}_i + \min\left(M_{i+1/2}, 0\right)\begin{bmatrix}
\rho a \\
\rho au \\
\rho av \\
\rho aw \\
\rho aH
\end{bmatrix}_{i+1}

M_{i+1/2} = M^+_i+M^-_{i+1}

M^\pm = \begin{cases}
\pm\frac{1}{4}\left(M\pm 1\right)^2 & (\left|M\right| \le 1) \\
\frac{1}{2}\left(M\pm\left|M\right|\right) & (\left|M\right| > 1)
\end{cases}

Der Druck wird wie folgt aufgeteilt.

p_{i+1/2} = p^+_i+p^-_{i+1}

p^\pm = \begin{cases}
\frac{p}{4}\left(M\pm 1\right)^2\left(2\mp M\right) & (\left|M\right| \le 1) \\
\frac{p}{2}\left(M\pm\left|M\right|\right)/M & (\left|M\right| > 1)
\end{cases}

5. Berechnungsbeispiel

5-1 Stoßwellenröhre

Lassen Sie uns als erstes Beispiel eine numerische Berechnung einer eindimensionalen Stoßwellenröhre durchführen. Dies ist die erste numerische Berechnung von komprimierbarem Fluid. Die Anfangsbedingungen sind wie folgt angegeben.

\rho_0 = \begin{cases}
1.29\ \mathrm{kg/m^3} & (x \ge 0\ \mathrm{cm}) \\
12.9\ \mathrm{kg/m^3} & (x < 0\ \mathrm{cm})
\end{cases} \\
u_0 = 0\ \mathrm{m/s}, T_0 = 300\ \mathrm{K}, \gamma = 1.4

Stellen Sie den Berechnungsbereich wie folgt ein.

0\ \mathrm{ms} \le t \le 5\ \mathrm{ms}, -50\ \mathrm{cm} \le x \le 50\ \mathrm{cm}

\Delta t = 0.005\ \mathrm{ms}, \Delta x = 0.5\ \mathrm{cm}

An der linken und rechten Grenze ist die Geschwindigkeit Null und der Temperatur- / Druckgradient Null. Die folgende Abbildung zeigt das Ergebnis der Berechnung mit der ungefähren Lehman-Lösung von AUSM und Roe.

5-2. Hindernisse im Luftstrom

Wenn ein Überschallstrom mit einem Objekt kollidiert oder wenn sich ein Objekt mit Überschallgeschwindigkeit bewegt, wird eine Stoßwelle um dieses herum erzeugt. Die Stoßwellen, die sich um Kampfflugzeuge und Kugeln bilden, sind berühmt. Als nächstes simulieren wir dies. Die Anfangsbedingungen sind wie folgt angegeben. Diesmal ist es zweidimensional.

\rho_0 = 1.29\ \mathrm{kg/m^3}, u_0 = v_0 = 0\ \mathrm{m/s}, T_0 = 300\ \mathrm{K}, \gamma = 1.4

Stellen Sie den Berechnungsbereich wie folgt ein.

0\ \mathrm{ms} \le t \le 100\ \mathrm{ms}, -100\ \mathrm{cm} \le x, y \le 100\ \mathrm{cm}

\Delta t = 0.01\ \mathrm{ms}, \Delta x = \Delta y = 1\ \mathrm{cm}

Platzieren Sie außerdem ein festes Objekt im Bereich von $ -5 \ \ mathrm {cm} \ le x, y \ le 5 \ \ mathrm {cm} $. Setzen Sie die Geschwindigkeit um das Objekt und an der Grenze von $ y = \ pm 100 \ \ mathrm {cm} $ auf Null und den Temperatur- / Druckgradienten auf Null. Als nächstes legt die Grenze von $ x = \ pm 100 \ \ mathrm {cm} $ die Randbedingung unter Verwendung der in 3-3 beschriebenen Lehman-Invariante fest. Insbesondere verwenden $ I_0, I_ + und I_- $, die sich von außerhalb der Grenze ausbreiten, voreingestellte Werte, und diejenigen, die sich von innen ausbreiten, werden unter Verwendung der Aufwinddifferenz berechnet. Diesmal ist $ u = 1,8a_0 = 624,94 \ \ mathrm {m / s}, \ rho = \ rho_0, T = T_0 $ als Werte außerhalb der Grenze von $ x = -100 \ \ mathrm {cm} $, $ x Geben Sie $ u = 0 \ \ mathrm {m / s}, \ rho = \ rho_0, T = T_0 $ als Werte außerhalb der Grenzen von = 100 \ \ mathrm {cm} $ an. Mit anderen Worten, der Luftstrom von Mach 1.8 strömt von der linken Seite. Die folgende Abbildung zeigt das Ergebnis der Berechnung mit AUSM. Sicherlich bilden sich Stoßwellen um das Objekt.

Lassen Sie uns abschließend prüfen, ob diese Stoßwelle korrekt berechnet wurde. Es ist bekannt, dass der folgende relationale Ausdruck ** Rankin-Yugonio ** vor und nach einer stationären vertikalen Stoßwelle gilt (Referenz [^ 8]).

\begin{eqnarray}
\frac{u_2}{u_1} &=& \frac{\rho_1}{\rho_2} = 
\frac{\left(\gamma+1\right)+\left(\gamma-1\right)p_2/p_1}{\left(\gamma-1\right)+\left(\gamma+1\right)p_2/p_1}\\
\frac{T_2}{T_1} &=& \frac{a_2^2}{a_1^2} = 
\frac{p_2}{p_1}\frac{\left(\gamma+1\right)+\left(\gamma-1\right)p_2/p_1}{\left(\gamma-1\right)+\left(\gamma+1\right)p_2/p_1}
\end{eqnarray}

Die Indizes 1 und 2 stehen stromaufwärts bzw. stromabwärts der Stoßwelle. Lassen Sie uns den durch die Berechnung erhaltenen Wert in diesen einfügen.

\frac{u_2}{u_1} = 0.41,\ 
\frac{\rho_1}{\rho_2} = 0.42,\ 
\frac{\left(\gamma+1\right)+\left(\gamma-1\right)p_2/p_1}{\left(\gamma-1\right)+\left(\gamma+1\right)p_2/p_1} = 0.42

\frac{T_2}{T_1} = 1.54,\ 
\frac{a_2^2}{a_1^2} = 1.54,\ 
\frac{p_2}{p_1}\frac{\left(\gamma+1\right)+\left(\gamma-1\right)p_2/p_1}{\left(\gamma-1\right)+\left(\gamma+1\right)p_2/p_1} = 1.53

Es ist ein gutes Gefühl!

[^ 8]: Tomomasa Tatsumi (1982) "Fluiddynamik" New Physics Series 21 Bakufukan ISBN 978-4-563-02421-5 (4-563-02421-X)