[PYTHON] Numerical calculation of compressible fluid by finite volume method
1.First of all
For flows around fast-moving objects such as aircraft and rockets, and for phenomena with strong pressure, temperature, and density changes such as combustion and explosion, it is necessary to consider the compressibility of the fluid. Compressible fluids cause various phenomena such as sound waves / shock waves and aerodynamic heating that are not found in incompressible fluids, and their numerical calculations have different difficulties than incompressible fluids. In this article, I would like to introduce a method for numerically calculating a compressible fluid using the Finite Volume Method (FVM).
2. Formulation of finite volume method
Consider the following partial differential equations for the vector $ \ mathbf {U} $, which is an array of $ N $ independent variables.
\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}
At this time, $ \ mathbf {F}, \ mathbf {G}, \ mathbf {H} $ is called ** flux **, and $ \ mathbf {S} $ represents the source term from the outside. Consider solving this equation numerically by the finite volume method. In the finite volume method, the calculation area is divided into cells called ** control volume **. This control volume does not have to be a uniform rectangular parallelepiped, and an unstructured grid can be used. As with the finite element method, one of the great strengths of the finite volume method is that it can handle complex shapes that match real problems. Now, let the above equation be the $ n $ th time $ t_n $ to the $ n + 1 $ th time $ t_ {n + 1} $ (assuming the time step is $ \ Delta t $) Integrate with the $ i $ th control volume $ V_i $ for space (reference [^ 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}
The time integral of one item on the left side can be performed immediately, and using the Gaussian divergence theorem, the volume integral of two items on the left side is converted to the surface integral on the surface $ \ partial V_i $ of the control volume.
\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}
However, $ \ mathbf {U} ^ n = \ mathbf {U} \ left (t_n, \ mathbf {x} \ right) $. If this surface $ \ partial V_i $ is composed of several more faces $ \ Omega_ {i, l} $ (for example, 6 rectangles for a rectangular parallelepiped), write the surface integral in the form of the sum for each face. I will.
\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}
Therefore,
\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}
By defining, we get:
\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
At this time, $ \ mathbf {F} ^ n_ {i, l}, \ \ mathbf {G} ^ n_ {i, l}, \ \ mathbf {H} ^ n_ {i, l} $ is ** numerical flow It is called a bundle **. This is the most common form of the finite volume method. If the control volume is evenly divided into the rectangular parallelepiped $ V_ {i, j, k} $ on each side of $ \ Delta x, \ Delta y, \ Delta z $, it can be transformed as follows.
\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
It has become quite familiar.
The subscripts $ \ mathbf {F} ^ n_ {i + 1/2, j, k} $ here are $ V_ {i, j, k} $ and $ V_ {i + 1, j, k} $ It means that the value of the boundary is used.
All of the above are mathematically exact variants, but in reality it is not possible to determine each flux on the left side exactly. Therefore, the essence of the finite volume method comes down to the question of how to approximate the numerical flux on this interface.
Below is a schematic diagram of the case where each control volume is rectangular in two dimensions. The value in the control volume is determined by the "balance" of the amount that goes in and out through the boundary surface.
3. Properties of compressible fluid
3-1. System of equations for compressible fluids
Now, let's sort out the system of equations for compressible fluids. By considering the conservation laws of density, momentum, and energy, the equation can be written in the form written at the beginning of 2 as follows (detailed derivation will be given to the textbook of fluid mechanics).
\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}
Here, $ \ rho and p $ are densities and pressures, $ u, v and w $ are velocities in the $ x, y and z $ directions, respectively, and $ E and H $ are energy and enthalpy per unit mass. Represents. $ \ Gamma $ is the specific heat ratio. Also, the viscosity / heat conduction and source terms are ignored here.
3-2. Diagonalization of Jacobian
An important difference between compressible and incompressible fluids is the presence of pressure waves or sound waves. Therefore, in addition to the information transmitted along the flow, the information is also transmitted by sound waves. So for each flux, the Jacobian $ \ mathbf {A} = \ partial \ mathbf {F} / \ partial \ mathbf {U}, \ mathbf {B} = \ partial \ mathbf {G} / \ partial \ mathbf Try to find {U}, \ mathbf {C} = \ partial \ mathbf {H} / \ partial \ mathbf {U} $. For example, for $ \ mathbf {A} $, it is calculated as follows (reference [^ 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}
Where $ a = \ sqrt {\ gamma p / \ rho} $, which represents the speed of sound. All of this Jacobian eigenvalues are real and can be diagonalized using the matrix $ \ mathbf {K} $ as follows (this calculation is very tedious, I don't want to do it again).
\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. Characteristic velocity and Lehman invariant
For simplicity, let's consider the one-dimensional case. The system of equations for one dimension is as follows.
\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}
At this time, by introducing the following amounts $ I_0, I_ \ pm $, we can show that the following equation holds.
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
Furthermore, assuming that $ I_0 $ is always a constant, that is, the ** isentropic condition **, $ I_ \ pm $ is simplified as follows.
I_\pm = u \pm \frac{2a}{\gamma-1}
This is called ** Lehman invariant **, which means that $ I_0, I_ +, and I_- $ propagate at velocities $ u, u + a, and u-a $, respectively. In other words, it was found that the true nature of the eigenvalues obtained in 3-2 was the propagation velocity of these invariants. Therefore, this eigenvalue is called ** characteristic velocity ** because it is the velocity that characterizes the system of equations.
In the case of one dimension, there are three independent variables, so if $ I_0, I_ +, I_- $ are found, all the variables can be found. In other words, it can be said that the value of each point is determined by the combination of the three types of waves that have propagated. The figure below schematically shows this.
4. Application to finite volume method
4-1. Upwind difference
To keep things simple, consider the following one-dimensional scalar advection equation.
\frac{\partial u}{\partial t} + \frac{\partial f}{\partial x} = 0,\ f = cu
Discretizing this by the finite volume method, it becomes as follows.
\frac{u_i^{n+1}-u_i^n}{\Delta t} + \frac{f_{i+1/2}^n-f_{i-1/2}^n}{\Delta x} = 0
At this time, $ c $ represents the speed of advection, and the method of assigning the flux value according to the sign of this is called the ** upwind difference method **.
f_{i+1/2}^n = \begin{cases}
cu_{i}^n & (c \ge 0) \\
cu_{i+1}^n & (c < 0)
\end{cases}
In the advection equation, the value is always determined by the information from the upstream, so the value is selected in this way. Now, how can we extend this for multivariable cases? As mentioned earlier, in the case of multivariables, there is a mixture of multiple types of waves, which must be well separated. Because it is possible that some waves propagate forward and some waves backward, it is necessary to change the direction of the upwind difference depending on the wave. Therefore, using the diagonalization considered in 3-2, we introduce a new vector $ \ mathbf {W} = \ mathbf {K} ^ {-1} \ mathbf {U} $. Assuming that the change in $ \ mathbf {K} $ is small, the equation can be transformed as follows.
\frac{\partial\mathbf{W}}{\partial t} + \mathbf{\Lambda}\frac{\partial\mathbf{W}}{\partial x} = \mathbf{0}
So 5 real eigenvalues $ \ lambda_1 = ua, \ lambda_2 = \ lambda_3 = \ lambda_4 = u, \ lambda_5 = u + a $ and $ \ mathbf {W} = \ left [W_1, W_2, W_3, W_4, W_5 \ right] $ transforms the equation into five independent transfer equations:
\frac{\partial W_l}{\partial t} + \lambda_l \frac{\partial W_l}{\partial x} = 0
Then, it seems good to apply the upwind difference to each of these five advection equations separately. So, if you introduce the amount $ \ lambda_l ^ + = \ max \ left (\ lambda_l, 0 \ right), \ lambda_l ^-= \ min \ left (\ lambda_l, 0 \ right) $, this advection equation will windward. The equations discretized by the above difference can be combined into one equation as follows.
\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
Then $ \ 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) $ and put this in the matrix Write 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}
Finally, $ \ mathbf {A} ^ + = \ mathbf {K} \ mathbf {\ Lambda} ^ + \ mathbf {K} ^ {-1}, \ mathbf {A} ^-= \ mathbf {K} \ As mathbf {\ Lambda} ^-\ mathbf {K} ^ {-1} $, the formula for $ \ mathbf {U} $ is as follows (using the assumption that the Jacobian change is small). The subscript of the position is deceived)
\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. Flux difference separation method
$\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,\ If you write \ Delta {\ mathbf {F} ^ +} ^ n_i + \ Delta {\ mathbf {F} ^-} ^ n_i = \ Delta \ mathbf {F} ^ n_i $, the last expression in 4-1 will be
\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}
It can be seen that the flux difference $ \ Delta \ mathbf {F} ^ n_i $ is separated in the positive and negative directions. The calculation method based on this formula is called ** Flux Difference Splitting (FDS) **.
Let's transform this a little more.
\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}
So, the 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}
Can be written as. Well, I've been fooling the Jacobian position subscripts all the time, but how do you actually calculate them? So Roe introduced a special mean weighted by the square root of the density as follows [^ 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}
With them
\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)
This is called ** Roe's approximate Lehman solution **.
4-3. Flux separation method
As another method, the flux is directly separated without going through the Jacobian as follows, and the Jacobian $ \ partial \ mathbf {F} ^ + / \ partial \ mathbf {U}, \ partial \ mathbf {F} ^ It is conceivable to configure the eigenvalues of-/ \ partial \ mathbf {U} $ to be non-negative and non-positive, respectively.
\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
This type of method is called ** Flux Vector Splitting (FVS) **.
4-3-1. Steger-Warming scheme
The simplest method is to separate the flux as follows, referring to the formula in 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
This is called the ** Steger-Warming scheme ** [^ 5]. In addition to this, various methods such as van Leer's flux separation method [^ 6] have been proposed.
4-3-2. AUSM Flux in the $ x $ direction
\mathbf{F} = \begin{bmatrix}
\rho u \\
\rho u^2+p \\
\rho uv \\
\rho uw \\
\rho uH
\end{bmatrix}
If you look at, you can see that everything except pressure is in the form of advection proportional to the flow velocity $ u $. Therefore, Liou and Steffen proposed a method called ** Advance Upstream Splitting Method (AUSM) **, which divides the flux into advection components and pressure components and separates them as follows [^ 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}
Here, $ M = u / a $ is called ** Mach number ** and is an important parameter when considering compressible fluids. In AUSM, the absolute value of this Mach number $ M $ is set to switch the expression at the boundary of 1 (that is, whether it is supersonic or not). First, for the advection component $ \ mathbf {F} ^ {\ left (a \ right)} $, do the following (the subscript n will be omitted for a while).
{\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}
The pressure is divided as follows.
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. Calculation example
5-1. Shock tube
As a first example, let's perform a numerical calculation of a one-dimensional shock tube. This is the first numerical calculation of compressible fluid. The initial conditions are given as follows.
\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
Set the calculation area as follows.
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}
At the left and right boundaries, the velocity is zero and the temperature / pressure gradient is zero. The figure below shows the result of calculating this with the approximate Lehman solution of AUSM and Roe.
Show video ·density  ·pressure  ·temperature  ·speed  ・ Mach number 
5-2. Obstacles in the airflow
When a supersonic stream collides with an object, or when an object moves at supersonic speed, a shock wave is generated around it. The shock waves formed around fighters and ammunition are famous. Next, let's simulate this. The initial conditions are given as follows. This time it is two-dimensional.
\rho_0 = 1.29\ \mathrm{kg/m^3}, u_0 = v_0 = 0\ \mathrm{m/s}, T_0 = 300\ \mathrm{K}, \gamma = 1.4
Set the calculation area as follows.
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}
In addition, place a fixed object in the range of $ -5 \ \ mathrm {cm} \ le x, y \ le 5 \ \ mathrm {cm} $. Set the velocity to zero and the temperature / pressure gradient to zero around the object and at the boundary of $ y = \ pm 100 \ \ mathrm {cm} $. Next, the boundary of $ x = \ pm 100 \ \ mathrm {cm} $ imposes the boundary condition using the Lehman invariant described in 3-3. More specifically, $ I_0, I_ +, and I_- $ that propagate from outside the boundary use preset values, and those that propagate from inside are calculated using the windward difference. This time, $ u = 1.8a_0 = 624.94 \ \ mathrm {m / s}, \ rho = \ rho_0, T = T_0 $ as values outside the boundary of $ x = -100 \ \ mathrm {cm} $, $ x Give $ u = 0 \ \ mathrm {m / s}, \ rho = \ rho_0, T = T_0 $ as values outside the boundaries of = 100 \ \ mathrm {cm} $. In other words, the airflow of Mach 1.8 flows from the left side. The figure below shows the result of calculating this using AUSM. Certainly, a shock wave is formed around the object.
Show video ·density  ·pressure  ·temperature  ・ Velocity in x direction  ・ Velocity in y direction  ・ Mach number 
Finally, let's check if this shock wave is calculated correctly. It is known that the following ** Rankine-Hugonieo relational expression ** holds before and after a stationary vertical shock wave (reference [^ 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}
Subscripts 1 and 2 represent the upstream and downstream of the shock wave, respectively. Let's put the value obtained by calculation into this.
\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
it is a good feeling!
[^ 8]: Tomomasa Tatsumi (1982) "Fluid Dynamics" New Physics Series 21 Baifukan ISBN 978-4-563-02421-5 (4-563-02421-X)
Recommended Posts