[PYTHON] Calcul numérique du fluide compressible par la méthode des volumes finis
1.Tout d'abord
Pour les écoulements autour d'objets en mouvement rapide tels que les avions et les fusées, et pour les phénomènes avec de forts changements de pression, de température et de densité tels que la combustion et l'explosion, il est nécessaire de considérer la compressibilité du fluide. Les fluides compressibles provoquent divers phénomènes tels que les ondes sonores, les ondes de choc et l'échauffement aérodynamique qui ne se trouvent pas dans les fluides incompressibles, et leurs calculs numériques présentent des difficultés différentes de celles des fluides incompressibles. Dans cet article, je voudrais présenter une méthode de calcul numérique du fluide compressible à l'aide de la méthode des volumes finis (FVM).
2. Formulation de la méthode des volumes finis
Considérez ce qui suit comme une équation différentielle partielle pour le vecteur $ \ mathbf {U} $, qui est un tableau de variables indépendantes $ N $.
\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}
A ce moment, $ \ mathbf {F}, \ mathbf {G}, \ mathbf {H} $ est appelé ** flux **, et $ \ mathbf {S} $ représente le terme source de l'extérieur. Envisagez de résoudre cette équation numériquement par la méthode des volumes finis. La méthode des volumes finis divise la zone de calcul en cellules appelées ** volumes de contrôle **. Il n'est pas nécessaire que ce volume de contrôle soit un carré régulier et des grilles non structurées peuvent être utilisées. Comme pour la méthode des éléments finis, l'une des grandes forces de la méthode des volumes finis est qu'elle peut gérer des formes complexes qui correspondent à des problèmes réels. Maintenant, laissez l'équation ci-dessus être le $ n $ ème temps $ t_n $ au $ n + 1 $ ème temps $ t_ {n + 1} $ (en supposant que l'incrément de temps est $ \ Delta t $) par rapport au temps. Intégrez avec le volume de contrôle $ i $ th $ V_i $ pour l'espace (référence [^ 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}
L'intégration temporelle d'un élément sur le côté gauche peut être effectuée immédiatement, et en utilisant le théorème de divergence gaussienne, l'intégration en volume de deux éléments sur le côté gauche est convertie en zone sur la surface $ \ V_i $ partiel du volume de contrôle.
\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}
Cependant, $ \ mathbf {U} ^ n = \ mathbf {U} \ left (t_n, \ mathbf {x} \ right) $. Si cette surface $ \ partial V_i $ est composée de plusieurs autres faces $ \ Omega_ {i, l} $ (par exemple, 6 rectangles pour un carré), écrivez l'aire dans la somme de chaque face. Je vais.
\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}
Donc,
\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}
Si défini comme, nous obtenons:
\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
À ce stade, $ \ mathbf {F} ^ n_ {i, l}, \ \ mathbf {G} ^ n_ {i, l}, \ \ mathbf {H} ^ n_ {i, l} $ est ** un flux numérique Cela s'appelle un bundle **. Il s'agit de la forme la plus courante de la méthode des volumes finis. Si le volume de contrôle est également divisé en carrés $ V_ {i, j, k} $ de chaque côté de $ \ Delta x, \ Delta y, \ Delta z $, il peut être transformé comme suit.
\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
C'est devenu assez familier.
Les indices $ \ mathbf {F} ^ n_ {i + 1/2, j, k} $ ici sont $ V_ {i, j, k} $ et $ V_ {i + 1, j, k} $ Cela signifie que la valeur de la limite est utilisée.
Tous les éléments ci-dessus sont des variantes mathématiquement exactes, mais en réalité, il n'est pas possible de déterminer exactement chaque flux sur le côté gauche. Par conséquent, l'essence de la méthode des volumes finis se résume à la question de savoir comment approximer le flux numérique sur cette interface.
Vous trouverez ci-dessous un diagramme schématique du cas où chaque volume de contrôleur est rectangulaire en deux dimensions. La valeur du volume de contrôle est déterminée par le «solde» de la quantité qui entre et sort à travers la surface limite.
3. Propriétés du fluide compressible
3-1. Système d'équation du fluide compressible
Ici, organisons une fois le système d'équation du fluide compressible. En considérant chaque loi de conservation de la densité, de la quantité de mouvement et de l'énergie, l'équation peut être écrite sous la forme écrite au début de 2 comme suit (la dérivation détaillée est laissée au manuel de dynamique des fluides).
\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}
Ici, $ \ rho et p $ sont la densité et la pression, $ u, v et w $ sont des vitesses dans les directions $ x, y et z $, respectivement, et $ E et H $ sont l'énergie et l'enthalpie par unité de masse. Représente. $ \ Gamma $ est le rapport de chaleur spécifique. De plus, les termes viscosité / conduction thermique et source sont ignorés ici.
3-2. Diagonalisation jacobienne
Une différence importante entre les fluides compressibles et non compressibles est la présence d'ondes de pression ou d'ondes sonores. Par conséquent, en plus des informations transmises le long du flux, les informations sont également transmises par des ondes sonores. Donc pour chaque flux, le Jacobien $ \ mathbf {A} = \ partial \ mathbf {F} / \ partial \ mathbf {U}, \ mathbf {B} = \ partial \ mathbf {G} / \ partial \ mathbf Essayez de trouver {U}, \ mathbf {C} = \ partial \ mathbf {H} / \ partial \ mathbf {U} $. Par exemple, pour $ \ mathbf {A} $, il est calculé comme suit (référence [^ 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}
Où $ a = \ sqrt {\ gamma p / \ rho} $, qui représente la vitesse du son. Toutes ces valeurs propres jacobiennes sont réelles et peuvent être diagonalisées en utilisant la matrice suivante $ \ mathbf {K} $ (ce calcul est très fastidieux, je ne veux pas le refaire).
\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. Vitesse caractéristique et invariant de Lehman
Par souci de simplicité, considérons le cas unidimensionnel. Le système d'équations pour une dimension est le suivant.
\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}
A ce moment, en introduisant les montants suivants $ I_0, I_ \ pm $, nous pouvons montrer que l'équation suivante est vraie.
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
De plus, en supposant que $ I_0 $ est toujours constant, c'est-à-dire que la ** condition d'entropie égale **, $ I_ \ pm $ est simplifiée comme suit.
I_\pm = u \pm \frac{2a}{\gamma-1}
Ceci est appelé ** invariant de Lehman **, ce qui signifie que $ I_0, I_ + et I_- $ se propagent aux vitesses $ u, u + a et u-a $, respectivement. En d'autres termes, il a été constaté que la vraie nature des valeurs propres obtenues en 3-2 était la vitesse de propagation de ces invariants. Par conséquent, cette valeur propre est appelée ** vitesse caractéristique ** car c'est la vitesse qui caractérise le système d'équations.
Dans le cas d'une dimension, il y a trois variables indépendantes, donc si $ I_0, I_ +, I_- $ sont trouvés, toutes les variables peuvent être trouvées. En d'autres termes, on peut dire que la valeur de chaque point est déterminée par la combinaison des trois types d'ondes qui se sont propagées. La figure ci-dessous le montre schématiquement.
4. Application à la méthode des volumes finis
4-1. Différence au près
Pour simplifier les choses, considérez l'équation de transfert scalaire unidimensionnelle suivante.
\frac{\partial u}{\partial t} + \frac{\partial f}{\partial x} = 0,\ f = cu
La dispersion de ceci par la méthode des volumes finis donne ce qui suit.
\frac{u_i^{n+1}-u_i^n}{\Delta t} + \frac{f_{i+1/2}^n-f_{i-1/2}^n}{\Delta x} = 0
À ce stade, $ c $ représente la vitesse de transfert, et la méthode d'attribution de la valeur de flux en fonction du signe de cela s'appelle la ** méthode de différence au vent **.
f_{i+1/2}^n = \begin{cases}
cu_{i}^n & (c \ge 0) \\
cu_{i+1}^n & (c < 0)
\end{cases}
Dans l'équation de translocation, la valeur est toujours déterminée par les informations en amont, la valeur est donc sélectionnée de cette manière. Maintenant, comment pouvons-nous étendre cela à plusieurs variables? Comme déjà mentionné, dans le cas de plusieurs variables, il existe plusieurs types d'ondes, il faut donc bien les séparer. Puisqu'il est possible que certaines vagues se propagent vers l'avant et certaines vagues vers l'arrière, il est nécessaire de changer la direction de la différence de vent en fonction de la vague. Par conséquent, nous allons introduire un nouveau vecteur $ \ mathbf {W} = \ mathbf {K} ^ {-1} \ mathbf {U} $ en utilisant la diagonalisation considérée en 3-2. En supposant que le changement de $ \ mathbf {K} $ est petit, l'équation peut être transformée comme suit.
\frac{\partial\mathbf{W}}{\partial t} + \mathbf{\Lambda}\frac{\partial\mathbf{W}}{\partial x} = \mathbf{0}
Par conséquent, 5 valeurs uniques réelles $ \ lambda_1 = ua, \ lambda_2 = \ lambda_3 = \ lambda_4 = u, \ lambda_5 = u + a $ et $ \ mathbf {W} = \ left [W_1, W_2, W_3, W_4, W_5 \ right] $ transforme l'équation en cinq équations de translocation indépendantes:
\frac{\partial W_l}{\partial t} + \lambda_l \frac{\partial W_l}{\partial x} = 0
Ensuite, il semble bon d'appliquer la différence au vent séparément à chacune de ces cinq équations de translocation. Donc, si vous introduisez le montant $ \ lambda_l ^ + = \ max \ left (\ lambda_l, 0 \ right), \ lambda_l ^ - = \ min \ left (\ lambda_l, 0 \ right) $, cette équation de translocation s'enroulera. Les équations discrétisées par la différence ci-dessus peuvent être combinées en une seule équation comme suit.
\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
Alors $ \ 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) $ Écrivez en forme.
\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}
Enfin, $ \ mathbf {A} ^ + = \ mathbf {K} \ mathbf {\ Lambda} ^ + \ mathbf {K} ^ {-1}, \ mathbf {A} ^ - = \ mathbf {K} \ Comme mathbf {\ Lambda} ^ - \ mathbf {K} ^ {-1} $, la formule pour $ \ mathbf {U} $ est la suivante (en supposant que le changement en jacobien est petit). L'indice de la position est trompé)
\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. Méthode de séparation par différence de flux de flux
$\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,\ Si vous écrivez \ Delta {\ mathbf {F} ^ +} ^ n_i + \ Delta {\ mathbf {F} ^ -} ^ n_i = \ Delta \ mathbf {F} ^ n_i $, la dernière expression de 4-1 est
\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}
On peut voir que la différence de flux de flux $ \ Delta \ mathbf {F} ^ n_i $ est séparée dans les directions positive et négative. La méthode de calcul basée sur cette formule est appelée ** Division de la différence de flux (FDS) **.
Transformons cela un peu plus.
\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}
Alors, le flux
\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}
Peut être écrit comme. Au fait, j'ai tout le temps trompé l'indice de la position de Jacobian, mais comment devrais-je le calculer? Par conséquent, Roe a introduit une moyenne spéciale pondérée par la racine carrée de la densité comme suit [^ 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}
Avec eux
\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)
C'est ce qu'on appelle la ** solution de Lehman approximative de Roe **.
4-3. Méthode de séparation du flux de flux
Une autre méthode consiste à séparer le flux directement sans passer par le jacobien comme indiqué ci-dessous, et le jacobien $ \ partial \ mathbf {F} ^ + / \ partial \ mathbf {U}, \ partial \ mathbf {F} ^ Il est concevable de configurer les valeurs uniques de- / \ partial \ mathbf {U} $ pour qu'elles soient respectivement non négatives et non positives.
\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
Ce type de méthode est appelé ** Flux Vector Splitting (FVS) **.
4-3-1. Schéma de réchauffement des bœufs
La méthode la plus simple consiste à séparer le flux comme suit, en se référant à la formule en 4-1.
\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
C'est ce qu'on appelle le ** schéma Steger-Warming ** [^ 5]. En plus de cela, diverses méthodes telles que la méthode de séparation de flux de van Leer [^ 6] ont été proposées.
4-3-2. AUSM Flux dans la direction $ x $
\mathbf{F} = \begin{bmatrix}
\rho u \\
\rho u^2+p \\
\rho uv \\
\rho uw \\
\rho uH
\end{bmatrix}
Si vous regardez, vous pouvez voir que tout sauf la pression a une forme de transfert proportionnelle à la vitesse d'écoulement $ u $. Par conséquent, Liou et Steffen ont proposé une méthode appelée ** Advance Upstream Splitting Method (AUSM) **, qui sépare le flux en une composante de transfert et une composante de pression comme suit [^ 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}
Ici, $ M = u / a $ est appelé ** nombre de Mach ** et est un paramètre important lorsque l'on considère les fluides compressibles. Dans AUSM, la valeur absolue de ce nombre de Mach $ M $ est définie de sorte que l'expression soit commutée à la limite de 1 (c'est-à-dire qu'elle soit supersonique ou non). Tout d'abord, pour le composant de transfert $ \ mathbf {F} ^ {\ left (a \ right)} $, procédez comme suit (l'indice n sera omis pendant un moment).
{\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}
La pression est divisée comme suit.
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. Exemple de calcul
5-1. Tube à ondes de choc
Comme premier exemple, effectuons un calcul numérique d'un tube à onde de choc unidimensionnel. Il s'agit du premier calcul numérique du fluide compressible. Les conditions initiales sont données comme suit.
\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
Définissez la zone de calcul comme suit.
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}
Aux limites gauche et droite, la vitesse est nulle et le gradient température / pression est nul. La figure ci-dessous montre le résultat du calcul avec la solution de Lehman approximative d'AUSM et de Roe.
<détails> 
·pression
·Température
·la vitesse
・ Numéro de Mach
5-2. Obstacles au flux d'air
Lorsqu'un flux supersonique entre en collision avec un objet, ou lorsqu'un objet se déplace à une vitesse supersonique, une onde de choc est générée autour de lui. Les ondes de choc formées autour des avions de combat et des balles sont célèbres. Ensuite, simulons cela. Les conditions initiales sont données comme suit. Cette fois, c'est bidimensionnel.
\rho_0 = 1.29\ \mathrm{kg/m^3}, u_0 = v_0 = 0\ \mathrm{m/s}, T_0 = 300\ \mathrm{K}, \gamma = 1.4
Définissez la zone de calcul comme suit.
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}
De plus, placez un objet fixe dans la plage $ -5 \ \ mathrm {cm} \ le x, y \ le 5 \ \ mathrm {cm} $. Réglez la vitesse sur zéro et le gradient température / pression sur zéro autour de l'objet et à la limite de $ y = \ pm 100 \ \ mathrm {cm} $. Ensuite, la limite de $ x = \ pm 100 \ \ mathrm {cm} $ impose une condition aux limites en utilisant l'invariant de Lehman décrit en 3-3. Plus spécifiquement, $ I_0, I_ + et I_- $ qui se propagent de l'extérieur de la limite utilisent des valeurs prédéfinies, et ceux qui se propagent de l'intérieur sont calculés en utilisant la différence au vent. Cette fois, $ u = 1.8a_0 = 624.94 \ \ mathrm {m / s}, \ rho = \ rho_0, T = T_0 $ comme valeurs en dehors de la limite de $ x = -100 \ \ mathrm {cm} $, $ x Donnez $ u = 0 \ \ mathrm {m / s}, \ rho = \ rho_0, T = T_0 $ comme valeurs en dehors des limites de = 100 \ \ mathrm {cm} $. En d'autres termes, le flux d'air de Mach 1,8 s'écoulera du côté gauche. La figure ci-dessous montre le résultat du calcul à l'aide d'AUSM. Certes, des ondes de choc se forment autour de l'objet.
<détails> 
·pression
·Température
・ Vitesse direction X
・ Vitesse direction Y
・ Numéro de Mach
Enfin, vérifions si cette onde de choc est correctement calculée. On sait que la ** expression relationnelle Rankin-Yugonio ** suivante est valable avant et après une onde de choc verticale stationnaire (référence [^ 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}
Les indices 1 et 2 représentent respectivement l'amont et l'aval de l'onde de choc. Mettons ici la valeur obtenue par le calcul.
\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
C'est un sentiment agréable!
[^ 8]: Tomomasa Tatsumi (1982) "Fluid Dynamics" Nouvelle série de physique 21 Bakufukan ISBN 978-4-563-02421-5 (4-563-02421-X)
Recommended Posts