Sod's Shock-tube Model

The next system of equations we consider is the Euler equations for gas dynamics

 

$$\rho_t + (\rho u)_x$$ = 0,

$$(\rho u)_t + (p+(\rho u)^2/\rho)_x$$ = 0, (10)

$$e_t + ((e+p)(\rho u)/p)_x$$ = 0,

with initial conditions from Sod's shock tube problem [15]

$$(\rho, (\rho u), e) = \cases{ (1.0, 0.0, 2.5), \qquad & if $x\le 0.5$ , \cr (0.125,0.0, 0.25), \qquad & if $ x > 0.5$ ,\cr }$$ (11)

where $\rho$ represents the gas density, u represents the gas velocity, and e represents the total energy of the gas. The gas pressure p satisfies

$$p = (\gamma-1)\left(e - \frac12\rho v^2\right),$$

where $\gamma$ stands for the adiabatic gas index. For an ideal gas, $\gamma=1.4$.

This problem has been solved successfully with different moving mesh methods by Winkler, et al.[17] and Dorfi and Drury [4]. The solution starts with a hot high density gas in the region $0\leq x\leq 0.5$, and a cold low density gas in the region $0.5<x\leq 1$. The gas is initially at rest. At t=0 the diaphragm separating the two regions is removed, causing a shock wave to propagate into the low density medium and a rarefaction wave into the high density medium. These two flow regions are separated by a contact discontinuity. At x=0 and x=1 we impose reflecting boundary conditions. Analytic solutions can be obtained for the early phases of the evolution. This problem exhibits several interactions of nonlinear waves, shock reflection, shock merging, the interaction of a shock with a contact discontinuity, and the reflection of a rarefaction wave.

The initial conditions (11) are discontinuous and cannot be used to generate an initial equidistributing mesh. Thus, we replace the initial conditions (11) by a continuous function

 

$$\rho \frac {(1.0+0.125)}{2} + \frac{(1.0-0.125)}{2}\tanh(200(x-0.5)),$$ =

u = 0.0, (12)

$$\frac {(2.5+0.25)}{2} + \frac{(2.5-0.25)}2 \tanh(200(x-0.5)),$$ e =

and compute an equidistributing mesh by our initial mesh generator. There is no viscosity term in Eq. (10). Our moving mesh solver has difficulty. Therefore, we add a small artificial viscosity on the right side of each equation in (11),

$$\varepsilon\rho_{xx}, \quad \varepsilon(\rho u)_{xx}, \quad \varepsilon e_{xx},$$ (13)

where $\varepsilon= \delta L/N$, $\delta = 10^{-2}$, L is the length scale, and N is the number of nodes. Numerical results show that these techniques do not affect the accuracy of the original problems. However, they make moving mesh methods work better.

The spatial discretization we choose is simple central differencing. A long rarefaction wave develops during the integration. If the arclength monitor is used, more points will be distributed to the rarefaction wave, which is not necessary. Hence we use the modified curvature monitor (5) in our computations. Numerical results show that this increases the computation time by only 1/3 but greatly improves the accuracy at the rarefaction corner, contact discontinuity and shocks. This evolution of the equations contains three phases: in the first phase, the shock and contact discontinuity propagates to the right; at $t\simeq 0.28$ the shock reflects from the wall (x=1), and the shock propagates to the left, the contact continues propagating to the right; at $t\simeq0.40$ the shock interacts with the contact discontinuity, and thereafter the shock continues propagating to the left and the contact discontinuity remains almost static. These three phases can also be observed from the space-time plot in Fig. 9.

  

Figure 9: Grid-point motion for the shock-tube problem (11). Output for every other node at each time step.

We also plotted the numerical solutions at times t=0.0,0.1, 0.2, 0.3, 0.4, 0.5.

  

Figure 10: Moving mesh method for shock-tube problem (11). N=81. Plotted at t=0.0,0.1,0.2,0.3,0.4,0.5.

In Fig. 10, the shock is resolved very well. The reference solution is computed with 400 fixed grid points and second order Godunov methods. The shock computed by our moving mesh method is much sharper than the reference solution. Indeed, the minimum spacing near the shock is O(10^-5) which cannot be achieved by 400 fixed grid points. However, the contact discontinuity is not resolved very well. This is because our numerical scheme cannot capture the contact discontinuity, and the shock has more weight than the contact discontinuity during the mesh distribution so that more points concentrate on the region of the shock. If more points are used, the results are better. However, the under-resolution of the contact discontinuity is still a problem unless some artificial compression is performed. We should point out that this problem is not related to the continuous initial condition we have chosen or the small viscosity terms we have added. If we make the initial condition sharper and the viscosity term smaller, the contact discontinuity will still be under-resolved.

As shown in Fig. 9, the time steps become very small at $t\simeq 0.28$, where the shock reflects from the wall (x=1), and $t\simeq0.40$, where the shock interacts with the contact discontinuity. This is because at that time, one discontinuity disappears suddenly and the grid points are pulled out to resolve the remaining structures; after that time, the discontinuity emerges and the grid points are pulled in again. We observed that from t=0.0 to t=0.2, it takes 72 time steps while from t=0.2 to t=0.3 it takes 105 steps. This irregularity of grid motion can be alleviated by the asymmetric time-filtering technique proposed by Winker, et al.[17]. This technique, however, is hard to incorporate into our moving mesh solver for general problems.