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 u*0= 10^{3}*m/s in the x-direction that separates the tangential velocities*
*v**L*= −5 × 10^{3}m/s on the left and*v**R*= 10^{3}m/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 p*_{0} = 10^{5} Pa in a two-dimensional shock tube of size [0*, 1] × [0, 1/10] m*^{2}and 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)}*= ρv*^{2}*/2, consistently by the method. Because of this, the*
pressure computed via (37) would not yield the accurate result that is in equilibrium with
*p*_{0} = 10^{5}Pa. 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,∂**t**q* *+ 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 energy*K** ^{(t)}*,

*∂**t**K*^{(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)*

with*K** ^{(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 solution

*K*

^{(t)}*6= (ρv)*

^{2}

*/(2ρ). Numerical*results present in Fig. 15 indicate that with the correction of

*K*

*to (42) this is a good*

^{(t)}*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,∂**t**q+ 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,*∂**t**q**+ 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 velocities*v**L* = 10^{3}m/s on the left and*v**R* = −5 × 10^{3}m/s
on the right of the interface. The computational domain is again a two-dimensional shock
tube of size [0*, 1] × [0, 1/10] m*^{2}*. 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 energy*K = ρ(u*^{2}*+ v*^{2}*)/2 from the jumps of the*
Riemann problem solution across the waves. Note that like the update of the state variables
*Q*^{n}* _{i j}*to the solutions of (36), in our wave propagation method, there is no problem to compute
the new value of

*K*

*i j*

*independently over a time step*

^{n}*1t 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 Q*

^{n}*and*

_{i j}*K*

^{n}*i j*is done, we may compute the pressure from the further revised formula of (42),

*p*^{n}_{i 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 of*K 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*, u*0*, −v*0*, p*0*), (ρ, u, v, p)*2*= (2ρ*0*, u*0*, v*0*, p*0*),*
*(ρ, u, v, p)*3*= (ρ*0*, −u*0*, v*0*, p*0*), (ρ, u, v, p)*4*= (3ρ*0*, −u*0*, −v*0*, p*0*),*

where *ρ*0= 10^{3} kg/m^{3}*, u*0= 10^{3} m/s, *v*0= 7 × 10^{2} *m/s, and p*0= 1 GPa. In this
prob-lem, the fluid in the first and third quadrants is a gas with*γ = 1.4, a = 1 Pa m*^{6}/kg, and
*b*= 10^{−4}m^{3}/kg, while the fluid in the second and fourth quadrants is a liquid with*γ = 4.4*
and*B = 6×10*^{8}Pa; 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, 3p*0*),* *(ρ, u, v, p)*2*= (ρ*0*, u*0*, ¯u, p*0*),*
*(ρ, u, v, p)*3*= (0.8ρ*0*, ¯u, ¯u, p*0*), (ρ, u, v, p)*4 *= (ρ*0*, ¯u, v*0*, p*0*),*

where*ρ*0= 1 kg/m^{3}*, u*0*= v*0*= 0, p*0= 10^{5} Pa, and ¯*ρ = 2.1 kg/m*^{3}*, ¯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*^{−3}m^{3}/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.