Theory¶

Finite Volume Methods for Euler’s Equation¶

In the present chapter, we briefly derive the discretized finite volume scheme for the invscid compressible Euler’s equations. The formulation is cell-centered with first order accuracy. The implementation of the same using OpenFOAM libraries is discussed. Finally, the parallelization of the solver using OpenFOAM library is discussed.

Governing Equations¶

Shown in equations (1)-(3) are the compressible Euler equations in conservative form. The momentum equation (2) is arranged such that the flow variables are on the left hand side and the driving potential terms are on the right.

(1) $\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0$

(2) $\frac{\partial (\rho \mathbf{u})}{\partial t}+ \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) = -\nabla p + \rho \mathbf{f}$

(3) $\frac{\partial (\rho e)}{\partial t}+ \nabla \cdot (\rho \mathbf{u} e ) = -\nabla \cdot (p \mathbf{u})$

The governing equations used for the CFD calculations are the compressible Euler’s equations in Cartesian coordinates. The equations were introduced in differential form in equations (1)- (3), but for the finite volume discretization, the integral form of the equations is employed. The integral equations in vector form is shown in equation (4).

(4) $\int_V{\dot{Q} dV} + \int_S{\mathbf{n}\cdot \mathbf{F} } = 0$

where,

(5) $Q = \left[ \begin{array}{c} \rho \\ \rho u \\ \rho v \\ \rho w \\ e \end{array} \right], \hspace{15pt} F_n = \mathbf{n \cdot F} = \left[ \begin{array}{c} \rho \bar{u} \\ \rho u \bar{u} +p n_x\\ \rho v \bar{u} +p n_y \\ \rho w \bar{u} +p n_z \\ e \bar{u} + p u_n \end{array} \right]$

The vector $\mathbf{n}(n_x,n_y,n_z)$ is the unit outward pointing normal of the surface of integration $S$ . The dot on top of any variable denotes differentiation with respect to time. The fluid velocity vector is $\mathbf{u}(u,v,w)$ and the grid velocity is denoted by $\mathbf{v}$ . The quantity $u_n=\mathbf{n \cdot u}$ and $v_n=\mathbf{n \cdot v}$ are the velocity components in the normal direction. Similarly, we have $\bar{u}=\mathbf{n \cdot (u - v)}=u_n-v_n$ .

Finite Volume Discretization¶

The volume integrals in equation (4) are evaluated using the cell center value $Q_c$ (piecewise constant) and the cell volume $V_c$ as shown in equation (6).

(6) $\int_V{QdV} \approx Q_c V_c$

Piecewise-constant reconstruction

In this piecewise-constant reconstruction one has to solve a local Riemann problem at every edge indicated by the fractional indices $i+1/2$, $i-1/2$ shown in figure Piecewise-constant reconstruction. The face average flux $F_n$ (piecewise-constant) is approximated as the face center value and it is obtained by solving the Riemann problem at the cell edge (given the left and right states). Introducing the above approximations, the discretized equation can be written as shown in equation :eq:euler_src_integral_discrete.

(7) $\dot{Q_c} + \sum_{i}^{n}{\frac{F_n S_i}{V_i}} = 0$

The flux $F_n$ can be evaluated using an Approximate Riemann solver for example the Roe scheme, using the left and right states of a face $Q_L$ and $Q_R$ as shown in equation (:eq:roe_scheme).

(8) $\hat{F}_{approx}(Q_L,Q_R,S)=\frac{1}{2}\left[\hat{F}(Q_L,S)+\hat{F}(Q_R,S) - \|\tilde{A}(Q_L,Q_R,S)\|(Q_R-Q_L)\right]$

The discrete equation shown in equation (7) can be rearranged and written in a semi-discrete form as an Ordinary Differential Equation (ODE) shown in equation (9). The operator $L$ consists of the discrete flux operator and the contribution due to source term. The explicit Euler time integration method has been implemented to solve this ODE as shown in equation (10).

(9) $\dot{Q} = L[Q]$

(10) $Q^{n+1} = Q^n + \frac{\Delta t}{v_{cell}} L[Q^n]$

where, $\Delta t$ is the time step and $v_{cell}$ is the cell volume. $L[Q^n]$ is the right hand side of the discrete equation (residue).

Local Time Stepping (LTS)¶

The solution is accelerated to steady state using the local time stepping (LTS). In equation (10) the $\frac{\Delta t}{v_{cell}}$ term is replaced by the CFL number ( $N_c$ ) and local maximum eigenvalue integral as shown in equation (11).

(11) $\frac{\Delta t}{v_{cell}} = \frac{N_c}{ \displaystyle\sum_{i=1}^{N_f} \left( |U_n| + a \right) dS }$

where, $N_f$ is the total number of faces in a cell and $dS$ is the face area of the $i^{th}$ face.

Boundary Condition¶

The following types of boundary conditions have been implemented, namely,

Supersonic Inflow

As the name suggests, these boundary conditions are imposed when the Mach number of the the flow entering or leaving the boundary of the computational domain is fully supersonic. Therefore, there is no influence of the downwind disturbances. One can make use of this property and prescribe the inflow/outflow boundary conditions. For a supersonic inflow boundary face, the flow parameters at the inlet are used to obtain the fluxes. For a supersonic outflow boundary face, the solution state extrapolated from within the interior computational domain is used for calculating the fluxes across that face.

Slip-wall

The slip wall boundary condition physically imposes a zero mass flux crossing the rigid wall. This can be written mathematically as $u_n = 0$ . This immediately suggests the following numerical flux formulation for the wall boundary face,

(12) $F_n^{wall} = \left[ \begin{array}{c} \rho u_n \\ \rho u u_n + p n_x \\ \rho v u_n + p n_y \\ \rho w u_n + p n_z \\ (e + p ) u_n \end{array} \right] = \left[ \begin{array}{c} 0 \\ p n_x \\ p n_y \\ p n_z \\ 0 \end{array} \right]$

The pressure value at the wall boundary face is usually computed by extrapolation from the interior computational domain. The nearest cell centroid value extrapolation has been implemented in the present work.

Riemann Extrapolation

The Riemann invariants for the 1D Euler’s equation (isentropic flow assumption) can be utilized as far-field boundary conditions by making a 1D flow assumption normal to the boundary face. If the far-field is sufficiently away from wall boundaries isentropic flow assumption is valid. The invariants $r_1,r_2$ and $r_3$ normal to a given finite volume interface is defined in equations shown below,

$r_1 = u_n + \frac{2 a }{\gamma - 1}, \hspace{5pt} r_2 = u_n + \frac{2 a }{\gamma - 1}, \hspace{5pt} and \hspace{5pt} r_3 = \frac{p}{ \rho^{\gamma} }$

One can detect the direction of propagation at the far-field boundary based on the sign of the Eigenvalues, which are defined as,

$\lambda_1 = u_n + a, \hspace{5pt} \lambda_2 = u_n - a, \hspace{5pt} and \hspace{5pt} \lambda_3 = u_n$

where, $u_n$ is the normal velocity to the interface and $a$ is the speed of sound at the interface. Based on the signs of the two variables four different cases can be considered as follows,

Supersonic Outflow ( $u_n \geq 0$ and $|u_n| \geq a$ )

The boundary condition is exactly same as that in previous section.

Supersonic Inflow ( $u_n < 0$ and $|u_n| \geq a$ )}

The boundary condition is exactly same as that in previous section.

Subsonic Outflow ( $u_n \geq 0$ and $|u_n| < a$ )

In this case, $\lambda_1 > 0, \hspace{5pt} \lambda_2 < 0 \hspace{5pt} and \hspace{5pt} \lambda_3 > 0$ , which means that two characteristic waves are outgoing and one wave is incoming. The boundary values are then determined as follows,

$\begin{eqnarray} u_n + a > 0 &\Longrightarrow& r_{1f} = r_{1e}\nonumber \\ u_n - a < 0 &\Longrightarrow& r_{2f} = r_{2 \infty }\\ u_n > 0 &\Longrightarrow& r_{3f} = r_{3e} \nonumber \\ u_{nf} &=& \frac{ r_{1e} - r_{2 \infty} }{ 2 } \\ a_f &=& \frac{\gamma - 1 }{4} ( r_{1e} - r_{2\infty}) \\ u_{\parallel f,1} &=& u_{\parallel e,1} \\ u_{\parallel f,2} &=& u_{\parallel e,2} \\ \rho_f &=& \rho_e \\ a_f &=& \sqrt{ \frac{\gamma p_f }{ \rho_f} } \end{eqnarray}$

Note that the subscripts $n, \hspace{5pt} \parallel, \hspace{5pt} f, \hspace{5pt} e$ and $\infty$ denote values normal/parallel to the boundary face, boundary cell values (unknowns), extrapolated value from interior and free-stream value respectively.

Subsonic Inflow ( $u_n < 0$ and $|u_n| < a$ )

In this case, $\lambda_1 > 0, \hspace{5pt} \lambda_2 < 0$ and $\lambda_3 < 0$ , which means that one characteristic wave is outgoing and two waves are incoming. The boundary values are then determined as follows,

$\begin{eqnarray} u_n + a > 0 &\Longrightarrow& r_{1f} = r_{1e}\nonumber \\ u_n - a < 0 &\Longrightarrow& r_{2f} = r_{2 \infty }\\ u_n < 0 &\Longrightarrow& r_{3f} = r_{3 \infty} \nonumber \\ u_{nf} &=& \frac{ r_{1e} - r_{2 \infty} }{ 2 } \\ a_f &=& \frac{\gamma - 1 }{4} ( r_{1e} - r_{2\infty}) \\ u_{\parallel f,1} &=& u_{\parallel e,1} \\ u_{\parallel f,2} &=& u_{\parallel e,2} \\ \rho_f &=& \rho_{\infty} \\ a_f &=& \sqrt{ \frac{\gamma p_f }{ \rho_f} } \end{eqnarray}$

Implicit non-reflecting boundary condition

Implicit boundary conditions are quite necessary for implementing Jacobian Free Newton-Krylov (JFNK) methods. Unlike the previous case, implicit non-reflecting boundary conditions fixes the normal boundary fluxes and not the boundary values. This way one ensures that the effect of boundary is implied in the residual calculation. Here one extrapolates the interior values for the left state and fix free-stream conditions as the right state and invokes the approximate Riemann flux solver to obtain the interfacial boundary fluxes and they are summed up to the neighboring cell residual. The implicit boundary described above is a weak boundary condition and hence is preferred more than the Riemann extrapolation due to the improved stability.

Navigation