Motivated by the work of Saurel and Abgrall [60] on slip lines (i.e., a contact discontinuity with jumps on the tangential velocities) for multicomponent problems, we now examine how our two-dimensional algorithm works for such a class of problems.
5.1. Preliminary
We begin by considering a simple numerical test in which the solution is a plan interface moving vertically by u0= 103m/s in the x-direction that separates the tangential velocities vL= −5 × 103m/s on the left andvR= 103m/s on the right of the interface. For convenience, we use the same two-phase flow setup as in Example 3.4.1 where the pressure is uniform with p0 = 105 Pa in a two-dimensional shock tube of size [0, 1] × [0, 1/10] m2and the fluid is a gas and liquid on the left and right of the interface, respectively. Note that this is a different kind of interface only problem that exists only in multiple space dimensions.
For this problem, we do the test using the first order wave propagation method with the Roe solver as usual; see Section 4.2. Results of a sample run with a 200×20 grid are shown in Fig. 14 at time t = 360 µs. From the figure, large errors in both the pressure and particle velocity in the x-direction are clearly seen. Carrying out more tests to the problem, we find that these errors remain at about the same order of magnitude as the mesh is refined and become even erroneous when a high-resolution method is employed instead. We note that this observation of the error is true also when we solve the problem using the shock-only Riemann solver. In fact, by following a similar analysis to that conducted in [60], we may explain the observed error behavior as being the failure to approximate the kinetic energy in the tangential direction,K(t)= ρv2/2, consistently by the method. Because of this, the pressure computed via (37) would not yield the accurate result that is in equilibrium with p0 = 105Pa. Note that this is a difficult problem to solve in practice and is even so in the single component case of the problem; see [80] for a general remark on numerical errors which occurred in this class of problem with Godunov-type schemes.
For this model slip line problem, analogously to the work done in [60], it is easy to improve upon the result by solving an extended one-dimensional system that combines the x-sweep of our two-dimensional model,∂tq + f (∂x, q) = 0, see (36), with an additional
FIG. 14. Preliminary results for a model slip line problem at time t= 360 µs; the first order method with the multicomponent model (36). The solid line is the exact solution and the points show the computed solution with a 200× 20 grid; only the solution along the cross-section of y = 1/20 m is presented.
transport equation for the tangential kinetic energyK(t),
∂tK(t)+ ∂x
¡K(t)u¢
= 0. (41)
Note that we have taken the above equation to be of the form that works even for shock and rarefaction waves as well. Rather than computing the pressure based on (37), which incurs large errors as shown in Fig. 14, here we use the following modified version,
p=
·
ρE −(ρu)2
2ρ − K(t)−
µγ − bρ γ − 1B
¶
−
µ2− γ − bρ γ − 1 aρ2
¶¸Áµ1− bρ γ − 1
¶ , (42)
withK(t)set by the solution of (41). Note that in the current slip line problem it is not dif-ficult to show the inconsistency of the numerical solutionK(t) 6= (ρv)2/(2ρ). Numerical results present in Fig. 15 indicate that with the correction ofK(t) to (42) this is a good approach without introducing any artificial oscillations in p and u, when performing the same test as in Fig. 14 with both the first order and high-resolution wave propagation methods. Observing the displayed profiles, the other components of the solutions are well behaved also.
It should be noted that no propagation of transverse waves and the solutions of the system in the y-direction,∂tq+ g(∂y, q) = 0, have been included in the above computation. This is eligible to do for this problem which is one-dimensional in nature. The Roe solver of this extended system can be derived quite easily as in Sections 3.2.2 and 4.2 and has been used to produce the results present in Fig. 15.
To further test the above corrected algorithm with the improvement of K(t) in the x-direction, we are next concerned with a two-phase liquid–gas Riemann problem studied by Saurel and Abgrall [60] in which the solution is composed of a leftward going rarefac-tion wave, a rightward going slip line, and a (weak) shock wave in front of the slip line. In this problem, we take the same initial data as concerned previously in Example 3.4.2, but
FIG. 15. Improved results of the test performed in Fig. 14; methods with the x-sweep of our two-dimensional model,∂tq+ f (∂x, q) = 0, and the additional equation for the tangential kinetic energy (41). Results using (a) the first order wave propagation method and (b) the high-resolution wave propagation method. The Roe solver is employed to the computations.
impose additionally tangential velocitiesvL = 103m/s on the left andvR = −5 × 103m/s on the right of the interface. The computational domain is again a two-dimensional shock tube of size [0, 1] × [0, 1/10] m2. Surely, the region in the y-direction is also redundant for this problem.
Figure 16 shows results for a run obtained using the high-resolution method with a 200× 20 grid up to time t = 240 µs, where the cross-section plots of ρ, u, v, and p along y= 1/20 m are presented at the final stopping time. From the displayed graphs, despite the slight overshoot inρ and undershoot in p at the tail of the rarefaction wave, see Fig. 4 also, we clearly observe the correct behavior of the computed slip line and also the rarefaction and shock waves as in comparison with the exact solution. We note that without using (42),
FIG. 16. High-resolution results for the Saurel–Abgrall liquid–gas Riemann problem at time t= 240 µs. The solid line is the exact solution and the points show the computed solution with a 200× 20 grid; only the solution along the cross-section of y= 1/20 m is presented.
the computed pressure obtained from (37) would become negative (which is nonphysical in the current case) within the first few time steps of the program execution. See [60] for the results of a similar calculation, where a MUSCL-type scheme is employed to solve the problem with the zero van der Waals gas constants: a= 0 and b = 0.
5.2. An Improved Multidimensional Algorithm
More generally, to be able to deal with multidimensional slip lines, shocks, and rarefaction waves at the same time, we should use the full set of the model system (36) and include suitable equations or schemes for the computation of the tangential kinetic energies in the x- and y-directions, respectively, yielding the accurate value for the pressure (cf. [60]).
Since the numerical method we have described in Section 4.2 is an unsplit one, after various unsuccessful attempts to generalize the approach proposed in Section 5.1, it turns out to be most convenient to implement the method based on the shock-only Riemann solver with an additional update to the total kinetic energyK = ρ(u2+ v2)/2 from the jumps of the Riemann problem solution across the waves. Note that like the update of the state variables Qni jto the solutions of (36), in our wave propagation method, there is no problem to compute the new value ofKi jn independently over a time step1t by the scheme. In fact, it is not difficult to show that the method reduces essentially to the scheme proposed in Section 5.1 for one-dimensional slip line problems, when the shock-only Riemann solver is employed there. When the update step of both Qni j andKni j is done, we may compute the pressure from the further revised formula of (42),
pni j+1=
·
ρE − K −
µγ − bρ γ − 1B
¶
−
µ2− γ − bρ γ − 1 aρ2
¶¸n+1 i j
Áµ1− bρ γ − 1
¶n+1 i j
, (43)
that is typically more accurate than simply employing (36) and (37) to the slip line problems (see comments and results shown below). Of course, in case there is not any strong shear flow moving along the interface, such as for a radially symmetric flow considered in Section 4.3, numerical results obtained using these two different approaches would be quite similar, where the solution ofK by the improved algorithm is approximately equal to [(ρu)2+ (ρv)2]/(2ρ) obtained from using the solutions of the basic conservation laws.
It should be mentioned that because the transfer of the energy between the kinetic and potential energies is governed implicitly by the conservation law of the total energy, in the general multidimensional case, we thus do not have a model equation for the motion of the total kinetic energy explicitly. For this reason, it is not clear at all how to form a linear system
of the governing equations and use the Roe solver as a basis to the current formulation of the scheme. Notably, this is a problem that we will work on in the future to develop an unsplit multicomponent algorithm for slip lines with a simpler approximate Riemann solver; this may not be an easy task, however. Two sample calculations are performed below to show the usefulness of the improved algorithm for slip line problems.
EXAMPLE5.2.1. As a first example, we consider a two-dimensional Reimann problem in which the initial condition is composed of four slip lines with the data in the four quadrants given by
(ρ, u, v, p)1= (ρ0, u0, −v0, p0), (ρ, u, v, p)2= (2ρ0, u0, v0, p0), (ρ, u, v, p)3= (ρ0, −u0, v0, p0), (ρ, u, v, p)4= (3ρ0, −u0, −v0, p0),
where ρ0= 103 kg/m3, u0= 103 m/s, v0= 7 × 102 m/s, and p0= 1 GPa. In this prob-lem, the fluid in the first and third quadrants is a gas withγ = 1.4, a = 1 Pa m6/kg, and b= 10−4m3/kg, while the fluid in the second and fourth quadrants is a liquid withγ = 4.4 andB = 6×108Pa; the domain is a unit square. We note that this problem is a multicompo-nent version of a case studied by Schulz-Rinne et al. where after breaking the membranes the slip lines spiral around its center in a clockwise manner forming an interesting vortex-like structure (Fig. 9 of [62]).
In Fig. 17, we present high-resolution results obtained using a 200× 200 grid and the improved method with the update of the total kinetic energy at time t= 150 µs. From the contours of the displayed quantities such as the density, particle velocities in the x- and y-direction, and pressure, we observe the very nice spiral structure of the computed so-lutions, without any artificial fluctuations near the slip lines. The cross-sectional plots of ρ, u, v, and p for the same run along the line x = 1 − y are shown in Fig. 18 where the solid line is the fine grid solution with a 400× 400 grid. We see good agreement between the two solutions.
FIG. 17. High resolution results for a two-dimensional Riemann problem, a slip-lines only case. The result is obtained using a 200× 200 grid and the improved method with the update of the total kinetic energy. Contours of the solutions for (a) density, (b) particle velocity in the x-direction, (c) particle velocity in the y-direction, and (d) pressure are shown at time t= 150 µs.
FIG. 18. Cross-sectional plots ofρ, u, v, and p for the same run shown in Fig. 17 along x = 1 − y. The solid line is the fine grid solution with a 400× 400 grid.
EXAMPLE5.2.2. Our second example is again a two-dimensional Riemann problem, but with a different set of data where the state in the first quadrant is connected to the second and fourth quadrants by a 1-shock moving leftward and downward of the x- and y-directions, respectively. The state in the third quadrant is connect to the second and fourth quadrants by a slip line, however. This problem is a variant of a run considered in Fig. 12 of [62]
where the interaction at the corner leads to a simple Mach reflection similar to what is seen in various shock reflections from boundaries. The data we use for the test are given by
(ρ, u, v, p)1= ( ¯ρ, ¯u, ¯u, 3p0), (ρ, u, v, p)2= (ρ0, u0, ¯u, p0), (ρ, u, v, p)3= (0.8ρ0, ¯u, ¯u, p0), (ρ, u, v, p)4 = (ρ0, ¯u, v0, p0),
whereρ0= 1 kg/m3, u0= v0= 0, p0= 105 Pa, and ¯ρ = 2.1 kg/m3, ¯u= 324.28 m/s for a Mach 1.65 shock wave. In this problem, with the exception that the fluid in the third quadrant is a polytropic gas withγ = 1.67, the fluid in the other quadrants is a constant covolume gas withγ = 1.4 and b = 10−3m3/kg.
We perform a similar test as in Example 5.2.1 using the improved high-resolution method, but with a finer 400× 400 grid for checking convergence of the solutions. Results of a sample calculation at time t= 800 µs are shown in Figs. 19 and 20. From the contours of the
FIG. 19. High resolution results for a two-dimensional Riemann problem, a case with simple Mach reflection.
The result is obtained using a 400× 400 grid and the improved method with the update of the total kinetic energy.
Contours of the solutions for (a) density, (b) particle velocity in the x-direction, (c) particle velocity in the y-direction, and (d) pressure are shown at time t= 800 µs.
FIG. 20. Cross-sectional plots ofρ, u, v, and p for the same run shown in Fig. 19 x = y.
plot, we observe sharp resolution of the primary shock waves and good behavior of the slip lines. The cross-sectional plots ofρ, u, v, and p for the same run along the line x = y provide an example of the basic structure of the solutions quantitatively. A detailed study of the algorithm to more general unstable interface problems such as the Kelvin–Helmholtz and Rayleigh–Taylor instabilities will be reported elsewhere in the future.
