_images/logo_small.png

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

_images/piecewise_const_reconst.png

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,

  1. 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.

  1. 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.

  1. 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}

  1. 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.