Grid-converged Solution and Analysis of the Unsteady Viscous Flow in a Two-dimensional Shock Tube

TL;DR

This paper provides a detailed, grid-converged numerical analysis of unsteady viscous flow in a 2D shock tube at Re=1000, revealing complex flow phenomena and establishing benchmark data using a high-order gas-kinetic scheme.

Contribution

It offers the first grid-converged solution for the shock tube problem at Re=1000 using a new high-order scheme, advancing numerical simulation accuracy for complex shock-driven flows.

Findings

Identification of vortex formation and shock bifurcation mechanisms

Detailed analysis of Kelvin-Helmholtz instability near contact surface

Benchmark data for high Reynolds number shock tube flows

Abstract

The flow in a shock tube is extremely complex with dynamic multi-scale structures of sharp fronts, flow separation, and vortices due to the interaction of the shock wave, the contact surface, and the boundary layer over the side wall of the tube. Prediction and understanding of the complex fluid dynamics is of theoretical and practical importance. It is also an extremely challenging problem for numerical simulation, especially at relatively high Reynolds numbers. Daru & Tenaud (Daru, V. & Tenaud, C. 2001 Evaluation of TVD high resolution schemes for unsteady viscous shocked flows. Computers & Fluids 30, 89-113) proposed a two-dimensional model problem as a numerical test case for high-resolution schemes to simulate the flow field in a square closed shock tube. Though many researchers have tried this problem using a variety of computational methods, there is not yet an agreed-upon grid-converged solution of the problem at the Reynolds number of 1000. This paper presents a rigorous grid-convergence study and the resulting grid-converged solutions for this problem by using a newly-developed, efficient, and high-order gas-kinetic scheme. Critical data extracted from the converged solutions are documented as benchmark data. The complex fluid dynamics of the flow at Re = 1000 are discussed and analysed in detail. Major phenomena revealed by the numerical computations include the downward concentration of the fluid through the curved shock, the formation of the vortices, the mechanism of the shock wave bifurcation, the structure of the jet along the bottom wall, and the Kelvin-Helmholtz instability near the contact surface.

Tables

Table 1. Table 1: Extracted data of the density along the bottom wall. \Rey = 200 \Rey 200 \Rey=200 .

$x$	$\rho$	$x$	$\rho$	$x$	$\rho$	$x$	$\rho$
0.3030	39.8418	0.6123	47.3367	0.7317	108.3916	0.8617	117.6452
0.4490	37.0662	0.6370	39.3203	0.7543	92.5760	0.9437	96.4287
0.5230	52.6465	0.6577	36.9558	0.7790	64.5319	0.9670	98.2689
0.5730	42.4400	0.6830	49.6513	0.7957	59.4386	0.9883	81.8465
0.5930	40.5506	0.7070	77.9810	0.8183	95.3607	0.9943	82.7077

Table 2. Table 2: Extracted data of the density along the bottom wall. \Rey = 1000 \Rey 1000 \Rey=1000 .

$x$	$\rho$	$x$	$\rho$	$x$	$\rho$	$x$	$\rho$
0.3001	39.5483	0.6063	36.6144	0.7491	72.8602	0.8817	118.8170
0.4391	36.9422	0.6173	37.3454	0.7579	56.5015	0.9081	106.8818
0.4525	38.1477	0.6405	30.9455	0.7621	67.0630	0.9239	102.2508
0.4811	60.1735	0.6581	35.1934	0.7701	74.9344	0.9447	97.2271
0.5085	53.2823	0.6703	45.5234	0.7839	53.1310	0.9631	99.4473
0.5121	53.4914	0.6761	42.6753	0.7909	50.8776	0.9739	77.4691
0.5265	44.3346	0.6891	67.0539	0.7991	65.2257	0.9785	81.3049
0.5355	43.1639	0.7111	94.2231	0.8051	85.1548	0.9829	77.2446
0.5631	38.6210	0.7391	69.4755	0.8121	82.7582	0.9923	97.4887
0.5769	38.9783	0.7451	67.8694	0.8399	106.8413	0.9999	85.6308

Equations

$$f_t+\mathbf{u}\cdot \nabla f=\frac{g-f}{\tau },$$

$$g=\rho {\left(\frac{\lambda }{\pi }\right)}^{\frac{K+2}{2}}{e}^{-\lambda \left[(u-U{)}^2+(v-V{)}^2+{\xi }^2\right]},$$

$$f(\mathbf{x},t,\mathbf{u},\xi )=\frac{1}{\tau }{\int }_0^{t}g\left({\mathbf{x}}^{\prime },t^{\prime },\mathbf{u},\xi \right)e^{-(t-t^{\prime })/\tau }d{t}^{\prime }+e^{-t/\tau }{f}_0\left(\mathbf{x}-\mathbf{u}t,\mathbf{u},\xi \right),$$

$$\tilde{g}(\mathbf{x},t,\mathbf{u},\xi )=g+g_{x}x+g_{y}y+g_{t}t+\frac{1}{2}{g}_{xx}{x}^2+\frac{1}{2}{g}_{yy}{y}^2+\frac{1}{2}{g}_{tt}{t}^2+g_{xy}xy+g_{xt}xt+g_{yt}yt.$$

$$f=g-\tau Dg=g-\tau \left(g_t+u{g}_x+v{g}_y\right).$$

$$\tilde{f}(\mathbf{x},t,\mathbf{u},\xi )=$$

$$+g_{xy}xy+g_{xt}xt+g_{yt}yt-\tau (g_t+g_{xt}x+g_{yt}y+g_{tt}t)$$

$$-\tau u(g_x+g_{xx}x+g_{xy}y+g_{xt}t)-\tau v(g_y+g_{xy}x+g_{yy}y+g_{yt}t).$$

$$f_0\left(\mathbf{x},\mathbf{u},\xi \right)=\{\begin{array}{rl}& f_0^l\left(\mathbf{x},\mathbf{u},\xi \right)={\tilde{f}}^l\left(g_0^l,\mathbf{x},0,\mathbf{u}\right),x\le 0,\\ & f_0^r\left(\mathbf{x},\mathbf{u},\xi \right)={\tilde{f}}^r\left(g_0^r,\mathbf{x},0,\mathbf{u}\right),x>0,\end{array}$$

$${\mathbf{W}}^l=\int g_0^l\mathbf{\psi }d\Xi ,{\mathbf{W}}^r=\int g_0^r\mathbf{\psi }d\Xi ,$$

$$\mathbf{\psi }={\left({\psi }_1,{\psi }_2,{\psi }_3,{\psi }_4\right)}^T={\left(1,u,v,\frac{u^2+v^2+{\xi }^2}{2}\right)}^T.$$

$$g\left(\mathbf{x},t,\mathbf{u},\xi \right)=\tilde{g}\left(g^e,\mathbf{x},t,\mathbf{u}\right),$$

$$\int g^e\mathbf{\psi }d\Xi ={\mathbf{W}}^e={\int }_{u\ge 0}{g}_0^l\mathbf{\psi }d\Xi +{\int }_{u<0}{g}_0^r\mathbf{\psi }d\Xi .$$

$$f(0,y,t,\mathbf{u},\xi )=$$

$$-e^{-t/{\tau }_n}\left[g^e-\left(u{g}_x^e+v{g}_y^e\right)t\right]+e^{-t/{\tau }_n}\left\{\begin{array}{rl}{g}^l& -\left(u{g}_x^l+v{g}_y^l\right)t,u\ge 0\\ g^r& -\left(u{g}_x^r+v{g}_y^r\right)t,u<0\end{array}\right\}.$$

$$\tau =\frac{\mu }{p^e},{\tau }_n=\tau +\alpha \Delta t{e}^{1-{\eta }^{-10}},\eta =\left|\frac{p^l-p^r}{p^l+p^r}\right|,$$

$$\mathbf{F}=\int uf\mathbf{\psi }d\Xi .$$

$${\mathbf{W}}_{ij}^{n+1}={\mathbf{W}}_{ij}^n$$

$$-\frac{1}{\Delta x_i\Delta y_j}{\int }_{t_n}^{t_{n+1}}{\int }_{-\frac{1}{2}\Delta x_i}^{\frac{1}{2}\Delta x_i}\left[{\mathbf{F}}_{j+1/2}(t,x)-{\mathbf{F}}_{j-1/2}(t,x)\right]dxdt.$$

$$a_x=g_x/g,a_y=g_y/g,a_t=g_t/g,$$

$$a_{xx}=g_{xx}/g,a_{yy}=g_{yy}/g,a_{xy}=g_{xy}/g,$$

$$a_{xt}=g_{xt}/g,a_{yt}=g_{yt}/g,a_{tt}=g_{tt}/g.$$

$$⟨\cdots ⟩=\int g(\cdots )\mathbf{\psi }d\Xi ,$$

$$⟨a_x⟩={\mathbf{W}}_x\to a_x,⟨a_y⟩={\mathbf{W}}_y\to a_y,⟨a_{x}u+a_{y}v+a_t⟩=0\to a_t,$$

$$⟨a_{xx}⟩={\mathbf{W}}_{xx}\to a_{xx},⟨a_{yy}⟩={\mathbf{W}}_{yy}\to a_{yy},⟨a_{xy}⟩={\mathbf{W}}_{xy}\to a_{xy},$$

$$⟨a_{xx}u+a_{xy}v+a_{xt}⟩=0\to a_{xt},⟨a_{xy}u+a_{yy}v+a_{yt}⟩=0\to a_{yt},$$

$$⟨a_{xt}u+a_{yt}v+a_{tt}⟩=0\to a_{tt}.$$

$$Q_i^r=\sum _{s=0}^2{w}_s{q}_s,Q_i^l=\sum _{s=0}^2{\tilde{w}}_s{\tilde{q}}_s,$$

$$q_0=\frac{1}{3}{\bar{Q}}_i+\frac{5}{6}{\bar{Q}}_{i+1}-\frac{1}{6}{\bar{Q}}_{i+2},{\tilde{q}}_0=\frac{11}{6}{\bar{Q}}_i-\frac{7}{6}{\bar{Q}}_{i+1}+\frac{1}{3}{\bar{Q}}_{i+2},$$

$$q_1=-\frac{1}{6}{\bar{Q}}_{i-1}+\frac{5}{6}{\bar{Q}}_i+\frac{1}{3}{\bar{Q}}_{i+1},{\tilde{q}}_1=\frac{1}{3}{\bar{Q}}_{i-1}+\frac{5}{6}{\bar{Q}}_i-\frac{1}{6}{\bar{Q}}_{i+1},$$

$$q_2=\frac{1}{3}{\bar{Q}}_{i-2}-\frac{7}{6}{\bar{Q}}_{i-1}+\frac{11}{6}{\bar{Q}}_i,{\tilde{q}}_2=-\frac{1}{6}{\bar{Q}}_{i-2}+\frac{5}{6}{\bar{Q}}_{i-1}+\frac{1}{3}{\bar{Q}}_i,$$

$$w_s=\frac{{\alpha }_s}{{\sum }_{p=0}^2{\alpha }_p},{\tilde{w}}_s=\frac{{\tilde{\alpha }}_s}{{\sum }_{p=0}^2{\tilde{\alpha }}_p},s=0,1,2,$$

$${\alpha }_s=\frac{d_s}{{\left(ϵ+{\beta }_s\right)}^2},{\tilde{\alpha }}_s=\frac{{\tilde{d}}_s}{{\left(ϵ+{\beta }_s\right)}^2},s=0,1,2,$$

$${\beta }_0=\frac{13}{12}{\left({\bar{Q}}_i-2{\bar{Q}}_{i+1}+{\bar{Q}}_{i+2}\right)}^2+\frac{1}{4}{\left(3{\bar{Q}}_i-4{\bar{Q}}_{i+1}+{\bar{Q}}_{i+2}\right)}^2,$$

$${\beta }_1=\frac{13}{12}{\left({\bar{Q}}_{i-1}-2{\bar{Q}}_i+{\bar{Q}}_{i+1}\right)}^2+\frac{1}{4}{\left({\bar{Q}}_{i-1}-{\bar{Q}}_{i+1}\right)}^2,$$

$${\beta }_2=\frac{13}{12}{\left({\bar{Q}}_{i-2}-2{\bar{Q}}_{i-1}+{\bar{Q}}_i\right)}^2+\frac{1}{4}{\left({\bar{Q}}_{i-2}-4{\bar{Q}}_{i-1}+3{\bar{Q}}_i\right)}^2,$$

$$d_0={\tilde{d}}_2=\frac{3}{10},d_1={\tilde{d}}_1=\frac{6}{10},d_2={\tilde{d}}_0=\frac{1}{10}.$$

$$\left(\rho ,u,v,p\right)=\{\begin{array}{rl}& (120,0,0,120/\gamma ),x\le 0.5,\\ & (1.2,0,0,1.2/\gamma ),x>0.5,\end{array}$$

$$ϵ=(f_2-f_1)/f_1.$$

Full text

1 Introduction

The shock tube is used as an experimental apparatus for studies of hypersonic flow and chemical reactions. The shock wave reflected from the end wall interacts with the boundary layer on the side wall induced by the incident shock as shown schematically in figure 1. Compression by the main high-energy flow from the left causes the fluid at the end wall to ‘leak’ backwards near the bottom wall where the fluid dynamic pressure is low because of the wall boundary layer. In time, the forward and backward flow in the boundary layer separates from the bottom wall resulting in a complex system of vortices, shock wave bifurcation, and other various flow structures. The homogeneity of the flow conditions in that region, however, is important for experimental tests using the shock tube (Bull & Edwards, 1968). Mark (1958) was the first to study this type of shock-wave/boundary-layer interaction. He developed a model based on the experimental results for analysis and prediction of the flow configuration and primary geometric parameters. Byron & Rott (1961) used a more realistic model, which is applicable for higher Mach numbers compared to Mark’s model. Subsequent theoretical analyses can be found in Davies & Wilson (1969) and Stalker & Crane (1978).

In recent decades, experiments and numerical simulations of this problem have been reported by other authors (Kleine et al., 1992; Wilson et al., 1995; Weber et al., 1995). As the viscosity plays an important role in the development of the flow field, the Reynolds number is a key parameter determining the features of the interaction. Differences of the Reynolds numbers used in the above papers make it difficult for comparison and analysis between their reported results.

Daru & Tenaud (2001) proposed a two-dimensional model problem for numerical simulation of the flow field in a viscous shock tube, which is designed for evaluating different numerical methods. This is a time-dependent unsteady problem. At moderate Reynolds numbers, a number of vortices appear in the computational domain due to high shearing effect, with length scales varying in a wide range. The multi-scale nature and the complicated flow field make it a good test case for high-order high-resolution schemes. As very fine grids are needed to resolve small structures, a practical problem is whether the computation could be completed within acceptable computational time. Therefore this case is a challenge for the robustness, accuracy, resolution as well as efficiency of a numerical method.

Since presented by Daru & Tenaud (2001), the viscous shock tube problem has been tested in many articles (Sjögreen & Yee, 2003; Daru & Tenaud, 2004; Kim & Kim, 2005a, b; Daru & Tenaud, 2009; Li et al., 2010; Houim & Kuo, 2011; Wan et al., 2012; Sun et al., 2014; Kotov et al., 2014; Tenaud et al., 2015; Wang & Ren, 2015; Pan & Xu, 2016; Pan et al., 2016). The cases with Reynolds numbers of 200 and 1000 are most frequently used. The results for the $\Rey=200$ case by different schemes are generally similar. But for the $\Rey=1000$ case, a range of solutions that are noticeably different have been reported in different papers. It seems that a grid-converged solution has not been shown at this Reynolds number. In this paper, grid-converged solutions are successfully obtained at both Reynolds numbers. The results for $\Rey=1000$ are in good agreement with the solution by Daru & Tenaud (2009). We recommend the current results be a reference solution.

The gas-kinetic scheme has been developed in the past years and shown great success in various categories of flows. The method employs the BGK equation (Bhatnagar et al., 1954) instead of the Navier–Stokes equations. A gas distribution function is modelled to represent the flow status. Since all macroscopic variables are simply the moments of the distribution function, the inviscid and viscous fluxes are treated simultaneously (Xu, 1998, 2001). Based on the high-order gas-kinetic scheme proposed by Luo & Xu (2013) which employs the WENO-JS reconstruction technique (Liu et al., 1994; Jiang & Shu, 1996) and a high-order gas evolution model, several simplifications are made by the authors and the resultant scheme enhances the efficiency by about $400\%$ for two-dimensional flows (Zhou et al., 2017). With this efficient high-order gas-kinetic scheme, we are able to simulate the viscous shock tube problem with finer grids and achieve grid-converged solutions at both $\Rey=200$ and $\Rey=1000$ in an acceptable CPU time.

In the following section, we will first outline the numerical method. §3 spells out the specification and computational conditions of the shock tube problem. The solutions at $\Rey=200$ and $\Rey=1000$ are presented in §4 and §5. §5 focuses on the difficult case at $\Rey=1000$. A procedure making use of the Grid-Convergence Index (GCI) is presented and used to prove grid convergence of our computations on a sequence of successively refined grids. The grid-converged solution provides fine details of the complex flow structure for the $\Rey=1000$ case. In §6 we discuss and analyse the detailed evolution of the fluid dynamics revealed by the numerical solution starting from the initiation of the incident shock wave and contact surface through a sequence of phenomena including the downward concentration of the fluid through the curved shock, the formation of the vortices, the bifurcation of the shock wave, creation of a jet-like flow towards the bottom wall, and vortex structures created by Kelvin-Helmholtz instability near the contact surface. Finally, we draw the conclusions in §7.

2 Numerical Procedure

In this section, we give a brief introduction to the numerical method. More details can be found in Luo & Xu (2013) and Zhou et al. (2017).

We start from the BGK equation (Bhatnagar et al., 1954):

[TABLE]

where $f$ is the gas distribution function, $g$ is the equilibrium state that $f$ approaches, $\bm{u}=(u,v)^{T}$ is the particle velocity, and $\tau$ is the collision time. For two-dimensional flow, the equilibrium (Maxwellian) distribution is

[TABLE]

where $\rho$ is the density, $U$, $V$ are macroscopic velocities in the $x$ and $y$ directions. $\lambda=m/2kT$, where $m$ is the molecular mass, $k$ is the Boltzmann constant and $T$ is the temperature. $K$ is the number of internal degrees of freedom which equals to $3$ for diatomic molecules. $\xi$ is the internal variable with $\xi^{2}=\xi_{1}^{2}+\xi_{2}^{2}+\cdots+\xi_{K}^{2}$.

(1) has an analytical integral solution:

[TABLE]

where $\bm{x}^{\prime}=\bm{x}-\bm{u}(t-t^{\prime})$ is the particle trajectory. Therefore $f$ depends on the equilibrium distribution function $g(\bm{x},t,\bm{u},\xi)$ and the initial distribution function $f_{0}(\bm{x},\bm{u},\xi)$.

Let $g=g(\bm{0},0,\bm{u},\xi)$ denote the Maxwellian distribution at the point $(x,y,t)=(0,0,0)$. Then $\tilde{g}$, the equilibrium distribution in the neighbourhood, can be expressed via the Taylor expansion to the second order:

[TABLE]

According to the Chapman-Enskog expansion, to the order of the Navier-Stokes equations, the non-equilibrium distribution $f$ has the following relation with the equilibrium distribution $g$ (Ohwada & Xu, 2004):

[TABLE]

Expand each term of $f$ at the point $(x,y,t)=(0,0,0)$, and neglect high-order derivatives of $g$, we have

[TABLE]

Note that for an arbitrarily given equilibrium state $g$, there exist $\tilde{g}$ and $\tilde{f}$ corresponding to $g$. Then we have the form $\tilde{g}=\tilde{g}(g,\bm{x},t,\bm{u})$, $\tilde{f}=\tilde{f}(g,\bm{x},t,\bm{u})$. The initial state at the cell interface should be discontinuous:

[TABLE]

where $g_{0}^{l}$ and $g_{0}^{r}$ correspond to the reconstructed conservative variables at the left and right sides of the cell interface, respectively, i.e.,

[TABLE]

where $d\Xi=dudvd\xi$, $d\xi=d\xi_{1}d\xi_{2}\cdots d\xi_{K}$, and $\bm{\psi}$ is the vector of moments:

[TABLE]

On the other hand, the equilibrium distribution function in the integral of the solution is replaced by

[TABLE]

where the equilibrium distribution $g^{e}$ is obtained from the statuses of both sides:

[TABLE]

Substitute the expressions of $f_{0}$ and $g$ into the solution (3), and neglect some unimportant terms (Zhou et al., 2017), the final form of the distribution function reads:

[TABLE]

The collision time is determined by

[TABLE]

where $\mu$ is the dynamic viscosity and $p^{e}$ is the pressure corresponding to $g^{e}$. $\tau_{n}$ is the numerical collision time which contains artificial dissipation (Luo & Xu, 2013). Note that an adaptive function $e^{1-\eta^{-10}}$ is designed for the numerical collision time. This function ensures that $\tau_{n}$ differs from $\tau$ only when the normalized pressure difference $\eta$ is large enough. By doing this we aim to provide a necessary but minimum artificial dissipation. $\alpha$ is a constant and is taken to be $0.3$ for all computations in this paper.

Once the distribution function $f$ is obtained, the flux at a vertically placed cell interface can be expressed as

[TABLE]

For a rectangular cell $[x_{i-1/2},x_{i+1/2}]\times[y_{j-1/2},y_{j+1/2}]$ with dimensions $\Delta x_{i}=x_{i+1/2}-x_{i-1/2}$ and $\Delta y_{j}=y_{j+1/2}-y_{j-1/2}$, the cell-averaged conservative variable $\bm{W}_{ij}$ is updated from the time $t_{n}$ to $t_{n+1}$ as follows:

[TABLE]

Since $\bm{F}$ is an explicit function of $t$ and $x,y$, the integrations in (15) can be easily obtained.

Finally, we give the coefficients for representing the derivatives of $g$ in (12):

[TABLE]

Each coefficient can be written as $\Lambda=\Lambda_{1}\psi_{1}+\Lambda_{2}\psi_{2}+\Lambda_{3}\psi_{3}+\Lambda_{4}\psi_{4}$. Define the moment of a variable as:

[TABLE]

then the coefficients are derived as follows:

[TABLE]

All moments can be obtained explicitly. See Xu (2001) for details.

To provide the initial values for the evolution process, the macroscopic variables and their derivatives need to be constructed before each computational step. In the perpendicular direction of the cell interface, a standard 5th-order WENO-JS method (Jiang & Shu, 1996) is used to determine the value of the variables on both sides of the interface. Following the suggestion in Shu (1997), the characteristic variables are used instead of conservative variables. For a scalar variable $Q$, assume $\bar{Q}_{i}$ is the averaged value in the $i$-th cell, $Q_{i}^{l}$ and $Q_{i}^{r}$ are the values to be reconstructed at the left and right boundaries of the $i$-th cell, then the process is as below:

[TABLE]

where

[TABLE]

We set $\epsilon=10^{-6}$ in our computations. The results of the one-dimensional WENO scheme are line-averaged values. A third-order interpolation is then used to obtain the value at the midpoint of the interface. After that, the first- and second-order derivatives in both $x$ and $y$ directions can be calculated from the reconstructed variables.

3 Description of the Viscous Shock Tube Problem

The viscous shock tube problem was proposed by Daru & Tenaud (2001). A diaphragm is vertically located in the middle of a square 2-D shock tube with unit side length, separating the space into the left and right parts. The initial state in non-dimensional form is given by

[TABLE]

where $\gamma=1.4$ is the specific heat ratio of air. The Prandtl number is taken to be $\Pran=0.73$. No-slip adiabatic conditions are applied at all boundaries of the tube.

The diaphragm is broken instantly at $t=0$. A shock wave with the Mach number $M\!a=2.37$ forms and moves towards the right, followed by a contact discontinuity. Simultaneously, a rarefaction wave expands in both directions. Figure 2 shows the evolution of density, velocity and pressure from $t=0$ to $t=1$ in the inviscid case (hence the flow is one-dimensional). It is seen from the figures that the incident shock reaches the right wall at about $t=0.21$. Then it is reflected back to the left, later interacting with the contact discontinuity.

With presence of viscosity, the incident shock wave induces boundary layers along the horizontal walls of the tube. They will then interact with the incident and reflected shock, as well as other structures appearing later. In figure 2, we can observe a number of wave reflections and interactions in the region close to the right wall for the one-dimensional inviscid case. For the two-dimensional viscous case, the flow field will surely be more complicated.

Since the configuration is symmetric about the line $y=0.5$, only half of the tube $[0,1]\times[0,0.5]$ is computed. And we focus on the evolution of the flow field from $t=0$ to $t=1$ at both Reynolds numbers of 200 and 1000. The viscosity is assumed to be constant (so that $\mu=1/\Rey$). All grids used are uniform with $\Delta x=\Delta y$. The CFL number is 1.0 for all computations.

4 The $\Rey=200$ Case

The $\Rey=200$ case has been simulated by many authors (Daru & Tenaud, 2001; Sjögreen & Yee, 2003; Daru & Tenaud, 2004; Kim & Kim, 2005a, b; Daru & Tenaud, 2009; Houim & Kuo, 2011; Wan et al., 2012; Sun et al., 2014; Kotov et al., 2014; Tenaud et al., 2015; Wang & Ren, 2015; Pan & Xu, 2016; Pan et al., 2016). At this relatively low Reynolds number, the results presented in different papers are quite consistent when the grid is fine enough. As reported in Daru & Tenaud (2009), the sufficient grid resolution is $1000\times 500$ for the high-order scheme presented therein. Other computations (Daru & Tenaud, 2001; Sjögreen & Yee, 2003) indicate the behaviour of high-order methods is obviously better than that of the second-order ones.

An important problem might be the lack of criteria for the judgement of convergence and for the comparison between results. Daru & Tenaud (2001, 2009) used the plot of density distribution along the bottom wall to demonstrate convergence. This method was also adopted by some other authors (Kim & Kim, 2005a, b; Pan et al., 2016). Another commonly used criterion is to compare the height of the primary vortex (Kim & Kim, 2005a, b; Wang & Ren, 2015; Pan & Xu, 2016; Pan et al., 2016). On the same uniform $500\times 250$ grid, the reported vortex height varies from 0.163 to 0.171 by different schemes. However, it is found that the flow structures are not necessarily the same even when the vortex heights are very close.

The grid convergence for the present scheme is illustrated in figure 3, where the density contours at $t=1$ are presented. The results by the $500\times 250$ grid and the $1000\times 500$ grid are almost indistinguishable. Figure 4(a) shows the density distribution along the bottom wall. The curves from the $500\times 250$ grid to the $1500\times 750$ grid are nearly identical. Even with a coarser $250\times 125$ grid, a very good result is obtained.

We think that the density distribution along the bottom wall is a good criterion for convergence study. Some critical points on the curve of the finest grid are extracted and listed in table 1 as a reference for comparison. The positions of the selected points are given in figure 4(b). For macroscopic evaluations of the computed results, we recommend the following three criteria which are easily measured in the density contour plot, see figure 5:

1. [(1)]

The position of the triple point, which is approximately $(x,y)=(0.58,0.137)$. 3. 2.

The height of the primary vortex, which is approximately 0.166. 4. 3.

The orientation of the long axis of the primary vortex. This is an obvious criterion for qualitative evaluation.

5 The $\Rey=1000$ Case: Numerical Simulation

The above case at $\Rey=200$ serves as verification for the present computational code. When the Reynolds number is increased to 1000, many fine flow structures appear hence the flow field becomes more complex. This case has also been simulated in several papers (Daru & Tenaud, 2001; Sjögreen & Yee, 2003; Daru & Tenaud, 2004, 2009; Li et al., 2010; Wan et al., 2012; Kotov et al., 2014; Pan et al., 2016). The results from different papers or even from different methods in the same paper are very different. One reason is the sensitivity of the problem to the computational conditions, another reason is that the grids used in previous studies are not fine enough to achieve grid convergence due to practical limit on computational time. Grid-convergence studies were performed in Daru & Tenaud (2001), Sjögreen & Yee (2003), and Daru & Tenaud (2009) with different numerical methods including classical TVD schemes and various high-order schemes. The most successful result is obtained by Daru & Tenaud (2009), where two high-order schemes (RK3-WENO5 and OSMP7) showed the same trend of convergence, and the results on the two finest grids ($3000\times 1500$ and $4000\times 2000$) are very similar. However, some small visible differences still exist on the two sets of grids, as noted in Daru & Tenaud (2009). Armed with the new accurate and efficient gas-kinetic scheme, we perform in this section a rigorous systematic grid-convergence study of the viscous shock tube problem at $\Rey=1000$.

5.1 Numerical results

Five successively refined grids are used for investigation, which are $1000\times 500$, $2000\times 1000$, $3000\times 1500$, $4000\times 2000$, and $5000\times 2500$, respectively. Figure 6 shows the density distribution at $t=1$ on different grids. It is clear that a converged solution in terms of the density field is obtained on the $3000\times 1500$ grid. And the main features of the vortex structures are able to be predicted on the $2000\times 1000$ grid.

The converged computational density distribution agrees well with the result on the finest $4000\times 2000$ grid in Daru & Tenaud (2009), providing evidence that the results obtained by Daru & Tenaud (2009) and by our current scheme are both accurate and reliable, thus can be regarded as a reference solution.

To perform a quantitative comparison, the density distribution along the bottom wall is shown in figure 7(a). The difference between the curves on the $2000\times 1000$ and $3000\times 1500$ grids is already very small. As in the $\Rey=200$ case, the critical points of the density distribution obtained on the finest grid are extracted and listed in table 2 as a reference. The positions of the selected points are shown in figure 7(b).

5.2 Grid refinement study with the Grid-Convergence Index approach

As shown in figure 6, we can hardly see any difference in the plot of the density distribution on the grid $3000\times 1500$, $4000\times 2000$, and $5000\times 2500$. Since the flow field is very complex, it is important to develop some quantitative measure on the convergence of the computational solutions to the presumed exact solution as the grid spacing is refined to approach zero. We adopt the Grid-Convergence Index (GCI) approach proposed by Roache (1994, 1997).

Based on the generalized theory of the Richardson Extrapolation (Richardson, 1911), the Grid-Convergence Index is defined to uniformly report the grid refinement tests. Assume $f_{1}$ and $f_{2}$ are solutions on a fine grid and a coarse grid, respectively, the relative error is expressed as

[TABLE]

Then the GCI of the fine-grid solution is defined by the following formula:

[TABLE]

where $r$ is the ratio of the grid spacing between the coarse and fine grids ($r=h_{2}/h_{1}>1$), and $p$ is the order of accuracy of the scheme. $c_{s}=3$ is a safety factor. As pointed out by Roache (1994), the GCI gives a conservative estimate of the error relative to the unknown ‘exact’ solution.

The underlying assumption of the GCI approach is the smoothness of the solution. The solution must have a Taylor series expansion at least up to the order of the numerical scheme. Despite the existence of many sharp ‘discontinuities’ in the present shock tube problem, the solution of the Navier-Stokes equations is not strictly discontinuous. Thus, the GCI still serves as a reliable measure on the convergence of our computations when the grids used are sufficiently fine enough.

In detail, the GCI on the $2000\times 1000$ and finer grids are computed. The calculations are performed on the target grid and the first coarser grid next to it, i.e., to get the GCI of the solution on the $3000\times 1500$ grid, the solutions on the $3000\times 1500$ grid and its neighbouring $2000\times 1000$ grid are used in (22) and (23).

In particular, we choose the $1000\times 500$ grid as a standard stencil. The GCI based on the averaged density in each stencil cell is computed. Since the cell numbers of all grids are integer multiples of the stencil cell number in both $x$ and $y$ directions, no interpolation or other approximation is needed. Following the suggestion of Roache (1994), since a uniform order $p$ can not be found all across the field which contains shocks and other discontinuities, a conservative value $p=1$ is used. After the GCI on each cell of the stencil is obtained, the average and root mean square of all the GCIs are taken and reported.

Roache (1994) also proposed a method for checking whether the asymptotic range of convergence is reached by two GCIs on three different grids, on the premise that the order of scheme is known. This is based on the fact that the GCI is essentially an estimate of the error level. Similarly, in the present case that the practical order of the scheme cannot be well defined, we assume that when three GCIs are located in a straight line in the log-log plot against the grid spacing, a conclusion can be drawn that the solution is converging with a constant order.

The results are shown in figure 8. For both averaging methods, the points corresponding to the $3000\times 1500$ grid, the $4000\times 2000$ grid and the $5000\times 2500$ grid are approximately in a line, indicating that the asymptotic range is achieved on the $3000\times 1500$ grid, whereas the result of the $2000\times 1000$ grid is out of the range. This conclusion agrees well with figure 6, where the visible details of the density distribution stay unchanged for the $3000\times 1500$ and finer grids, but not for the $2000\times 1000$ one.

If we go back to the original meaning of the GCI, it is seen in figure 8, from an overall perspective, that the averaged relative error of the result obtained by the $5000\times 2500$ grid is less than 1%, with respect to the exact solution.

The viscous shock tube problem at $\Rey=1000$ is naively simple in geometry and initial and boundary conditions. Yet, it encompasses the evolution of almost all elementary flow phenomena of a viscous compressible flow and their mutual interactions, resulting in a complex dynamic flow field with a multitude of fine scales. As such it offers a difficult but arguably necessary test case to demonstrate the accuracy and efficiency of modern high-resolution and high-order numerical methods for compressible viscous flows. The grid-converged solution for this problem as well as the rigorous GCI approach presented here provide the research community a useful database and approach in comparing and assessing different numerical methods for their numerical discretization, flux models, shock capture strategies, effect of numerical dissipation, time evolution, and implementation of boundary conditions.

6 The $\Rey=1000$ Case: Analysis of the Complex Flow Physics

The dynamic evolution process of the flow field at $\Rey=1000$ is of great significance for understanding the fluid dynamics of the interactions between boundary layers, vortices, and wave systems in supersonic flow. Analysis and discussion of the flow physics of this problem, however, has been rather minimum in previous papers except those by Daru & Tenaud (2004, 2009). (Chen (2015) calculated a slightly different problem and gave some discussions on the flow behaviours at early stages.) This is partly due to the complexity of this problem and partly lack of adequate proof of numerical convergence. With the solution on the $3000\times 1500$ grid proven above be grid converged, we proceed to present and analyse the details of the flow field and its time evolution. Important observations during the process are emphasised.

Before detailed description, we present the whole history of the physical dynamic process in figure 9, where the magnitude of the density gradient at different time points of interest are shown in chronological order.

At $t=0$, break of the diaphragm results in three different waves: a right-moving shock wave, a contact discontinuity following the shock, and an expansion wave propagating in both directions. The waves travel freely into the undisturbed region creating a boundary layer on the bottom wall behind. See figure 10. This configuration is similar to the inviscid case in figure 2, except for the creation of the wall boundary layer and thickening of the two discontinuities (especially the contact discontinuity) due to viscous effect.

The boundary layer is attached to and dragged by the right-moving shock wave, as can be seen in figure 10(b), where the distribution of the velocity in the $x$-direction is shown. The boundary layer thickens as one moves away from its initiation point at the foot of the shock much like a usual boundary layer over a flat plate until $x=0.75$ where the contact discontinuity is located. The effective Reynolds number is increased due to the high density in the freestream flow behind the contact surface, resulting in a decrease of the boundary layer thickness.

At this stage the boundary layer is behind the shock wave and is theoretically of zero thickness at the foot of the shock. Therefore, the shock front remains effectively straight across the channel and curves only slightly as it touches the wall. On the contrary, the contact discontinuity, being a material wave front that moves with the fluid, is dramatically bent over the boundary layer because of the no-slip condition on the wall. It is seen from figure 10(a) that a very oblique contact discontinuity is stretched along the horizontal wall and it connects with the vertical one outside the boundary layer.

The curved near-wall section of the shock wave gets enlarged with time. Since the pressure gradient is perpendicular to the shock surface, the curving of the shock generates a non-zero $y$-direction component of the pressure gradient. Figure 11(a) shows the distribution of the magnitude of the pressure gradient at $t=0.15$. The shock is more curved at locations closer to the wall. The $y$-component of the pressure gradient is shown in figure 11(b). Obviously this quantity is closely related to the curvature of the shock. As a consequence, the fluid will experience a sudden acceleration when it flows across the narrow shock, obtaining a downward velocity. Although this velocity is very small and nearly invisible because of the much larger flow velocity in the $x$-direction, later we will see that it is of great importance in the following dynamic process.

At about $t=0.21$, the right-travelling shock wave encounters the end wall and is then reflected by it. As the shock is curved, it reaches the wall successively from upper parts to lower parts. Figure 12 presents three snapshots around the time of reflection. In figure 12(a), the upper part shown in the plot has just moved to the wall; in figure 12(b), the upper part has been reflected back while the lower part just touches the wall; in figure 12(c), the lower part has also completed the reflection. Since theoretically the horizontal velocity of the flow in the region behind a reflected normal shock is zero, the downward-concentrating effect of the curved shock can be observed very obviously in figure 12(b) and 12(c). It is clear from the streamlines that the fluid flows to the lower-right corner from upper regions behind the reflected shock wave. However, we emphasise that this process started from the very beginning: A region with negative velocity in the $y$-direction always exists after the shock wave is generated, see figure 13(a). The gathering of flow near the root of the shock makes the density there larger, as shown in figure 13(b). To get a better view, a Galilean transform is made at $t=0.2$: A constant is subjected from the $U$ velocity in the flow field, so that the $V$ velocity is shown more clearly. The streamlines after transformation are presented in figure 14. It demonstrates how the fluid is moving to the bottom wall. This process has no essential difference with the phenomenon behind the reflected shock shown in figure 12.

We will then focus on the flow in the lower right corner. It is seen from figure 13(c) that the shock wave disperses near the bottom wall due to viscous effect. Hence it is more like a sequence of compressible waves in this region. In addition, the shock is very curved there and the strength in the $x$-direction is then weakened. As a consequence, the reflected wave in the near-wall region is not as strong as that in the upper region where the incident shock is thin and normal to the right wall. This effect creates a pressure gradient pointing to the lower left direction, see figure 15. Driven by such a pressure gradient, the downward flow alters its direction to the left. Figures 15(a) to 15(d) display the process how the streamlines adjust to the perpendicular direction to the pressure contour lines.

The reversed flow at the corner shown in figure 15(d) soon encounters the incident flow around the position of the left edge of the reflected shock. With continuous supply of fluid, an oblique separation line forms and gets longer between the two parts of the fluid. This process is shown in figure 16. In the last three snapshots of figure 16 we can see that the fluid beside the separation line is forced to flow downward or upward, generating two sink points at the ends of the separation line and a saddle point in the middle.

With the lifting of the upper sink point, its distance to the bottom wall increases, hence the fluid around the sink point has larger velocity and momentum. In this situation, the streamlines roll up forming a vortex around the point, which gets larger in size with entrainment of more fluids, see figure 17. Notice that the streamlines and the pressure contour lines finally adjust to be orthogonal to each other.

It is interesting that there is a close connection between the vortex and the oblique reflected shock wave. Notice that the left edge of the vortex is aligned with the oblique shock. The rotation of the vortex makes the difference on the left and right sides of the oblique shock larger so that the strength of the shock is enhanced. And the asymmetric pressure distribution in the direction parallel to the oblique shock caused by the vortex rotation makes the shock more oblique, as shown in figure 18. On the other hand, after the flow passes the oblique shock, the normal component of the velocity decreases to near zero, while the tangential component remains unchanged. Therefore, the fluid behind the oblique reflected shock flows upwards along it, which is right in the same direction with the rotating flow in the vortex. This means that the oblique shock provides a momentum injection mechanism to the vortex and makes it larger and stronger.

The process in this stage can also be interpreted in another view: The downward moving fluid behind the reflected shock wave carries higher momentum than the fluid in the boundary layer. Then it is easy for the former to insert inside the boundary layer, as shown in figure 19, where the momentum vectors and the distribution of the momentum magnitude are plotted.

At about $t=0.27$, the reflected shock wave encounters the right-travelling contact discontinuity and is nearly stopped by it. The contact discontinuity then moves on with a lower speed. Simultaneously, a new shock wave is formed and propagates to the right. The interaction process is presented in figure 20. Notice that the contact discontinuity has not reached the reflected shock wave in figure 20(a).

The flow in the bulk region outside the viscous boundary layer is similar to the one-dimensional inviscid case. The viscous flow in the near-wall region has a very different behaviour. Since the shock wave becomes oblique in the lower region, the status change of the flow passing the normal shock and the oblique shock is different. This difference of the two regions behind the reflected shock becomes extremely distinct after the contact discontinuity brings the large-density and high-momentum fluid behind it. Remember that the vortex is carrying fluid along the oblique shock from the lower region to the upper region. To accommodate the huge difference of the fluid property, a shock appears at the interface between the two regions, i.e., bifurcation occurs at the junction point of the normal shock and the oblique shock. See figure 20(c). This process is more clearly presented in figure 21, where we can see a lambda-shaped structure around the triple point.

At about $t=0.32$, the new shock wave produced by the shock/contact-discontinuity interaction has been reflected back by the right wall. It then crosses the right-moving contact discontinuity and is slowed down by it. After that, the shock interacts with the vortex and then with the stationary shock, making it start to move again to the left, along with the triple point of the lambda-shaped shock. There are also many other secondary waves and a number of interactions between them at this stage. But they are relatively weak hence do not affect the primary picture much.

Later when the vortex is stronger, it dominates the local flow field. We can see from figure 22 that the dense fluids are entrained by the vortical flow around the core of the vortex, creating a jet inserting into the bottom lighter fluids. The momentum magnitude distributions are plotted in figure 23, showing how the jet is generated at the lower right corner of the high-momentum region.

In another view, the jet is enclosed by two contact discontinuities, one of which originates from the vertical contact discontinuity while the other originates from the oblique contact discontinuity. This mechanism is clearly shown in figure 9(a), 9(b) and 9(c). In figure 9(a), the two contact discontinuities with different orientations are presented in the density-gradient-magnitude contour map. Then the horizontal discontinuity encounters the oblique shock wave and the vertical contact discontinuity encounters both the normal and oblique shocks. The two contact discontinuities become stronger after getting through the shock wave, and their shape remains the same, except that the horizontal one is a little deflected up by the oblique shock. Then they are both bent and carried down by the vortex, forming the two boundaries of the jet.

As the horizontal contact discontinuity is deflected behind the oblique shock wave, a wedge-shaped area appears between it and the bottom wall. In figure 24(a), we can see that the jet becomes longer and extends to the left, alternatively reflecting on the two boundaries of the wedge-shaped area. This area is then divided by the jet into several individual regions distributed on both sides of the jet. Small secondary vortices are induced by the jet in these individual regions. And these vortices may further induce smaller vortices, see the section between $x=0.85$ and $x=0.9$ in figure 24(a). This demonstrates the multi-scale feature of the flow field. To avoid ambiguity, the vortex formed at the beginning will be called the primary vortex hereinafter. It should be noted that a large vortex is generated by the primary vortex in the lower right corner.

From figure 24(a) which shows the distribution of the momentum magnitude, it is also found that the upper right edge of the primary vortex is a thin contact surface. Therefore the Kelvin-Helmholtz instability occurs around it, which is shown in figure 9(f), where a sequence of vortical structures is observed near the contact surface. These structures are driven by the primary vortex down to the corner, and merged with the stationary contact discontinuity located at around $x=0.94$. After getting to the bottom wall, the vortical structures are taken over by the anticlockwise-rotating corner vortex shown in figure 24(a). The corner vortex carries these structures upward along its streamlines. This process is presented in figure 9(g), 9(h) and 9(i). Meanwhile, the wide stationary contact discontinuity at about $x=0.97$ is rolled up. In figure 9(i), we can see that the big rotating structure at the lower right corner involves at least four contact discontinuities altogether.

The distribution of the momentum magnitude at $t=0.75$ is plotted in figure 24(b). It is clear that besides the left-moving zigzagging jet beneath the deflected contact discontinuity, there is another jet turning right at about $x=0.92$ along the streamlines of the corner vortex. In fact, this flow pattern can be found at each point where the jet impinges on the bottom wall, which is also the lower boundary of the wedge-shaped area: The jet is split by the wall into two branches due to its high momentum, the bigger turning to the left, and the smaller to the right (see figure 25).

As for the upper boundary of the wedge-shaped area, when the jet impinges on it, a part of the jet will be ejected up leaking into the outer region above the deflected contact discontinuity. This phenomenon is demonstrated by the streamlines in figure 24(b). Clearly that the jet is not totally restricted in the wedge-shaped area. The ejected fluids are then taken to the right by the high-momentum flow in the outer region, producing a pulling force which makes the jet broken at the contacting points, as is shown in figure 25, where a gap is seen at about $x=0.7$.

With presence of the gaps, the fluids under the jet are carried up by the ejected jet into the outer flow region. These hot and light fluids are also taken away by the outer high-momentum flow, forming thin filaments, the biggest among which finally bumps onto the left edge of the primary vortex. See the temperature distribution in figure 26.

On the other hand, the secondary vortices above the jet are lift up as the ejected part of the jet is taken to the right by the outer fluid. Under the flushing of the high-speed outer flow, they are deformed rapidly and get closer to the neighbouring vortices, shown in figures 27(a) to 27(d). Then the adjacent small vortices are merged into a big one since they share the same rotating direction, see figures 27(e) and 27(f). The rotation of the new big vortex tends to make the jet become straight. Also notice that the amount of the fluid under the jet has decreased due to the ejection from the gaps. The final result is that the small vortices in the wedge-shaped area all vanish, and the jet becomes very flattened.

7 Conclusion

The viscous shock tube problem is simulated by an efficient high-order gas-kinetic scheme. Grid-convergence studies by using the GCI approach are presented for the two cases at $\Rey=200$ and $\Rey=1000$. Grid-converged solutions are achieved on the $500\times 250$ grid for the $\Rey=200$ case and on the $3000\times 1500$ grid for the $\Rey=1000$ case. Nevertheless, critical points on the curve of the density distribution along the bottom wall are extracted from the result obtained on the finest grid ($1500\times 750$ for $\Rey=200$ and $5000\times 2500$ for $\Rey=1000$) as benchmark data. The viscous shock tube problem is a good test case for accuracy, resolution and efficiency of high-order high-resolution schemes. We hope the present results can serve as a benchmark solution.

The dynamic process of the $\Rey=1000$ case is analysed. Important evolution stages, flow structures and physical phenomena are interpreted in detail, including the downward flow due to the shock curvature, the origin of the primary vortex, the shock bifurcation (the formation of the lambda-shaped shock), the Kelvin-Helmholtz instability, the components of the corner vortex, the secondary vortices and their breaking up. These processes demonstrate the complexity of the interactions between shock waves, contact discontinuities, boundary layers, and multi-scale vortices.