Application of WENO-Positivity-Preserving Schemes to Highly Under-Expanded Jets

Abstract

The starting transient of highly under-expanded supersonic jets is studied by means of very high resolution weighted essentially non oscillatory finite volume schemes, coupled with a positivity-preserving scheme in order to ensure positivity of pressure and density for high compression/expansion ratio. Numerical behaviour of the schemes is investigated in terms of grid resolution, formal accuracy and different approximated Riemann solvers. The transient flow field is also discussed.

Appendices

Appendix

Code Verification and Validation

The numerical software used for the present work, the Open source Finite volume Fluid dynamics (OFF) code [28], is freely available at https://github.com/szaghi/OFF.

Five 1D numerical tests were considered for code verification and validation. In the first, a problem with a smooth solution was simulated, in order to assess the accuracy of the schemes used in the present paper. The other four are discontinuous problems: three are well known classical 1D Riemann’s problems, while the last one is a more discriminating benchmark, consisting of two interacting blast waves.

1.1 Problem with Smooth Solution

For the accuracy assessment of the 7th-order WENO scheme with PPS limiter used for the present paper, we solved the Euler equations with the initial conditions

$$\begin{aligned} \left\{ \begin{array}{l} \rho = 1 + \frac{4}{5}sin(\pi x) + \frac{1}{10}sin(5 \pi x) \\ u = \frac{1}{2}\left( x - \frac{1}{2}\right) ^4 \\ p = 10 + 2 x ^4 \end{array} \right. \end{aligned}$$

(15)

where the 1D domain (x) is discretized with a uniform grid with $N_i$ finite volumes. The specific heat ratio considered is $\gamma = \frac{c_p}{c_v}=1.4$. Differently from the under-expanded jet transient simulations above presented, for this test, a 10th order Runge–Kutta scheme with 17 stages [2]was used for time integration, in order to focus on space accuracy.

The initial conditions 15 are built to provide a smooth solution for the whole simulation time [0, 5], as shown in Fig. 16a.

Four uniform grid with increasing resolution were used, namely $N_i= [100, 200, 400, 800]$, the refinement ratio being kept constant and equal to 2. Figure 16b reports the density distribution the simulation time $t = 2.0$.

Fig. 16

Smooth solution test: density field. a Time history evolution in time interval [0,4.5], b grid convergence analysis at time 2.0

Given the $L_2$-norms of the differences between the density profiles on the three finest grids, the apparent convergence order is given by [17]

$$\begin{aligned} p = \frac{\textit{log}\left( ||\rho _{\textit{medium}}- \rho _{\textit{coarse}}||_2 / ||\rho _{\textit{fine}}- \rho _{\textit{medium}}||_2 \right) }{\textit{log}(r)} \end{aligned}$$

(16)

where $r=2$ is the refinement ratio. The computed value was $p=6.55$, the maximum difference between the solutions for the two finest grids being smaller than $10{-7}$.

1.2 Sod’s and Lax’s Problems

In all test cases, the x axis is discretized with a uniform grid of $N=100$ cells and the left and right boundary conditions are set as non-reflective.

The first test problem considered is the Sod’s problem [22], whose initial conditions are:

$$\begin{aligned} P_L = \left[ \begin{array}{l} \rho \\ u \\ p \end{array} \right] = \left[ \begin{array}{l} 1.000 \\ 0.000 \\ 1.000 \end{array} \right] \quad P_R = \left[ \begin{array}{l} \rho \\ u \\ p \end{array} \right] = \left[ \begin{array}{l} 0.125 \\ 0.000 \\ 0.100 \end{array} \right] \end{aligned}$$

(17)

where the specific heat ratio is $\gamma = c_p/c_v=1.4$. After the break of the initial discontinuity, a shock moving to the right and a rarefaction fan moving to left are produced, together with a contact discontinuity moving to the right. Figure 17 shows a comparison of the numerical solutions computed by means of the WENO-positivity-preserving scheme coupled with the Exact (iterative), the HLLC and LLF approximated solvers. The exact solution is also reported. The (non-dimensional) time unit considered is $t=0.2$. The four plots report the comparison for schemes with 1st, 3rd, 5th and 7th orders. For all accuracy orders the solution of the Exact solver and the HLLC one are very close, while the solution of the LLF solver is more dissipative (especially for low order approximations). Nevertheless, for high order schemes (${\ge } 5{th}$) the LLF solution is comparable with the Exact and HLLC ones. This is more clear when analyzing Fig. 18a, where the LLF solutions obtained with all schemes are compared with the exact solution. As the order of accuracy increases, the higher dissipation of LLF is mitigated by means of the enlarged reconstruction stencil. As Fig. 18b shows, the same behaviour is observed for HLLC solutions, but the differences between low and high order approximations are smaller. This test highlights the fact that the LLF solver, when coupled with high order scheme (${\ge } 5{th}$), can give accurate results, comparable with the solutions obtained with more expensive solvers (like the Exact and HLLC solvers). However, this is not true in general: for more discriminating tests, the dissipation of LLF solver introduces strong errors, especially when the flow field has strong discontinuities (e.g. in the two interacting blast waves test that follows). For multidimensional problems the inaccuracies due to LLF solver were found to be less relevant.

Fig. 17

Sod’s problem: comparison of exact, HLLC and LLF solvers with exact solution ($t=0.2$). a 1st order, b 3rd order, c 5th order, d 7th order

Fig. 18

Sod’s problem: solutions with different order of accuracy ($t=0.2$). a LLF solutions with different order of accuracy, b HLLC solutions with different order of accuracy

A more interesting test is the modified Sod’s problem [25]. A uniform translational velocity is added to the left state of the classical Sod’s problem so that the resultant rarefaction wave is transonic. This test is useful in assessing the entropy property of the scheme. The initial conditions are the following:

$$\begin{aligned} P_L = \left[ \begin{array}{l} \rho \\ u \\ p \end{array} \right] = \left[ \begin{array}{l} 1.000 \\ 0.750 \\ 1.000 \end{array} \right] \quad P_R = \left[ \begin{array}{l} \rho \\ u \\ p \end{array} \right] = \left[ \begin{array}{l} 0.125 \\ 0.000 \\ 0.100 \end{array} \right] \end{aligned}$$

(18)

where the specific heat ratio is $\gamma = \frac{c_p}{c_v}=1.4$.

Figure 19 reports a comparison of the density profiles for the exact, HLLC and LLF solvers with 1st and 7th order accuracy. The same plots for entropy are reported in Fig. 20. Exact and HLLC solvers have good accuracy even at first order, while the absence of contact discontinuity of the LLF introduces too much dissipation, and, when coupled with the first order scheme, it fails to capture the entropy profile. However, when coupled with high order reconstructions, it recovers the gap with respect to the other solvers and, at least for the 7th order scheme, the entropy profile is correctly captured.

Fig. 19

Modified Sod’s problem: density profiles with different solvers and order of accuracy ($t=0.2$). a 1st order, b 7th order

Fig. 20

Modified Sod’s problem: entropy profiles with different solvers and order of accuracy ($t=0.2$). a 1st order, b 7th order

Consider now the Lax’s problem [8]:

$$\begin{aligned} P_L = \left[ \begin{array}{l} \rho \\ u \\ p \end{array} \right] = \left[ \begin{array}{l} 0.445 \\ 0.698 \\ 3.528 \end{array} \right] \quad P_R = \left[ \begin{array}{l} \rho \\ u \\ p \end{array} \right] = \left[ \begin{array}{l} 0.500 \\ 0.000 \\ 0.571 \end{array} \right] \end{aligned}$$

(19)

where the specific heat ratio is again $\gamma = c_p/c_v=1.4$. The solution is made of a right traveling shock wave, a contact discontinuity and a left rarefaction fan. Figure 21 shows the density, pressure and velocity profiles of the numerical solutions computed by means of the HLLC solver for the all the four different accuracy orders. The (non-dimensional) time unit considered is $t=0.13$. As the accuracy increases, the resolution of the density profile wave is clearly improved.

Fig. 21

Lax’s Problem: HLLC solutions with different order of accuracy ($t=0.13$). a Density profiles, b velocity profiles

1.3 Interacting Blast Waves

The last 1D test case is the interaction of two blast waves [26]. The initial conditions are slightly different from the standard shock-tube problems above:

$$\begin{aligned}&P\left( {x,t = 0} \right) = {P_0}\left( x \right) = \left\{ \begin{array}{l} {P_L}\qquad 0.0< x< 0.1\\ {P_M}\qquad 0.1< x< 0.9\\ {P_R}\qquad 0.9< x < 1.0 \end{array} \right. \end{aligned}$$

(20)

$$\begin{aligned}&P_L = \left[ \begin{array}{l} \rho \\ u \\ p \end{array} \right] = \left[ \begin{array}{l} 1.00 \\ 0.00 \\ 1000.00 \end{array} \right] \, P_M = \left[ \begin{array}{l} \rho \\ u \\ p \end{array} \right] = \left[ \begin{array}{l} 1.00 \\ 0.00 \\ 0.01 \end{array} \right] \, P_R = \left[ \begin{array}{l} \rho \\ u \\ p \end{array} \right] = \left[ \begin{array}{l} 1.00 \\ 0.00 \\ 100.00 \end{array} \right] \qquad \quad \end{aligned}$$

(21)

where the specific heat ratio is, as usual, $\gamma = c_p/c_v=1.4$. For this test the space domain [0, 1] was discretized with uniform grid of $N=200$ finite volumes and reflective boundary conditions were used at both the left and right boundary. The (non-dimensional) time unit considered is $t=0.038$. This test is more discriminating than the previous ones, because the two strong blast waves develop and collide, producing a new contact discontinuity.

Figure 22 shows the computed solutions with Exact, HLLC and LLF solvers for all accuracy orders. The generated contact discontinuity at about $x=0.725$ is not well resolved for orders up to the 3rd. Higher order schemes (${\ge } 5{th}$) are able to resolve the contact discontinuity, but the LLF resolution is very poor if compared with the Exact and HLLC ones (Fig. 22d). Besides, the LLF solver dissipates too much the strength of the blast waves. This result suggests that, also when coupled with very high order schemes, the LLF solver can introduce too much dissipation.

Fig. 22

Interacting blast waves ($t=0.038$). a Exact solver solutions with different order of accuracy, b HLLC solutions with different order of accuracy, c LLF solutions with different order of accuracy, d Comparison of different Riemann’s problem solvers at orders 7th