Modified algorithm of the Riemann problem solution

1. Introduction

The use of front tracking methods for the modeling of a three-phase flow in a porous medium was initially difficult due to the fact that there are no analytical solutions to the corresponding Riemann problem, a particular case of the Cauchy problem, when the initial condition is two constants separated by a single break [1]. Most of the pre-existing solutions of the Riemann problem with respect to three-phase filtration were limited to the consideration of simplified functions of relative permeability, in which the relative permeability for each phase depends only on the saturation of this phase. But even when this restriction was weakened, and the relative permeabilities could be functions of all saturation, in published works on this topic, there were only general guidelines for constructing the solution (see, for example, [2]). Only relatively recently in [3-5] a complete set of solutions of such a problem was presented under the condition that certain limitations are imposed on the function of relative permeability. Moreover, the paper proposed effective algorithms for finding solutions. They were based on a correction-correction strategy in conjunction with the iterative Newton-Raphson procedure, which provided quadratic convergence. Since the solution of the Cauchy problem usually requires the solution of a large number of Riemann problems, the presence of a general and highly efficient solver of Riemann is topical. An example of such a solver is described in [1], which is also based on the use of the predictor-corrector strategy, the algorithm of which is as follows:

1. Set the left and right states: ${\mathbf{S}}_l$ and ${\mathbf{S}}_r$.

2. Set the initial approximation and the trial solution: ${\mathbf{S}}_m^{tr}$, $W_1^{\mathrm{t}\mathrm{r}}=R_1$, $W_2^{tr}=R_2$.

3. Find the solution for a given configuration and update the wave structure: $\left[{\mathbf{S}}_m,W_1,W_2\right]=\mathrm{W}\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{S}\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{c}\mathrm{t}\left({\mathbf{S}}_l,{\mathbf{S}}_r,{\mathbf{S}}_m^{tr},W_1^{tr},W_2^{tr}\right)$.

4. Verify the admissibility (monotonic increase in wave velocities): if ${\sigma }_1>{\sigma }_2$, then set a new initial approximation ${\mathbf{S}}_m$ and designate the solution as incorrect $W_1^{tr}=W_2^{tr}=$ 0.

5. Check the convergence of the algorithm: if $W_1{W}_2=W_1^{tr}{W}_2^{tr}$, then stop, otherwise assign $W_1^{tr}{W}_2^{tr}=W_1{W}_2$ and ${\mathbf{S}}_m^{tr}={\mathbf{S}}_m$. Go to step 3.

Given the left and right states ${\mathbf{S}}_l$ and ${\mathbf{S}}_r$, accordingly, the algorithm solves the Riemann problem by determining the intermediate state ${\mathbf{S}}_m$ and the types of waves $W_1$, $W_2$ connecting the three constants. The algorithm begins by setting the initial approximation ${\mathbf{S}}_m^{tr}$ of the intermediate state and assuming that the trial solution $R_1{R}_2$ is, i.e. consists of two rarefaction waves. This choice ensures that the expected intermediate state ${\mathbf{S}}_m$ will be inside the triangle of saturation [1]. The core of the algorithm is Section 3, in which the following two actions are performed:

1. An intermediate state is determined for a given trial wave structure $W_1^{tr}{W}_2^{tr}$ and initial approximation ${\mathbf{S}}_m^{tr}$.

2. It is found out which structure of the solution waves could be for the newly calculated intermediate state. The wave type is displayed separately for each of the families ($k=$ 1,2). Although in order to obtain a reliable implementation, this step must be developed with caution, the idea is quite simple: the admissible discontinuity of the $k$ family must satisfy the Lax entropy condition; if the discontinuity is not permissible, and the constant states connected with the help of the wave of the $k$ family are located on one side relatively $V_k$, then we have a rarefaction wave of the $k$ family, otherwise a mixed wave.

Due to the fact that for some saturation states the solution must include detached branches of the phase Hugoniot set, it is not sufficient to verify the admissibility of only individual waves. If there are gaps in the solution, then you need to make sure that they form an increasing sequence of wave velocities, i.e. that ${\sigma }_1<{\sigma }_2$ [1]. The algorithm finishes work if the trial wave structure is admissible $W_1^{tr}{W}_2^{tr}$. Otherwise, the values of the intermediate state and wave structure are updated to the calculated values. It is believed that the rarefaction and discontinuity curves in the phase space are usually located close enough to each other, so the position of the intermediate state is not very sensitive to the type of solution, and the procedure converges, as a rule, in one iteration [1].

2. Modeling results of the algorithm improvement for Riemann problem solving

Developed in [1, 3], the algorithms for solving the Riemann problem with respect to one-dimensional three-phase filtering in most cases yield fairly fast convergent results with a given accuracy. However, when testing the proposed approaches by setting randomly inside the triple diagram of left and right states, some limitations of the algorithm described above and possible ways of overcoming them were revealed.

We consider the case ${\mathbf{S}}_l={\left[\begin{array}{ll}\text{0,4}& \text{0,6}\end{array}\right]}^T$ and ${\mathbf{S}}_r={\left[\begin{array}{ll}\text{0,05}& \text{0,9}\end{array}\right]}^T$ for which in Fig. 1 shows the basic sets and curves that form the solution.

Fig. 1Integral curves and Hugoniot sets for Sl=0,40,6T and Sr=0,050,9T

It can be seen that the Hugoniot phase set for a state ${\mathbf{S}}_r$ includes a detached branch, and the integral curve and the local Hugoniot curve of the second family below the same state differ significantly from each other and actually are on opposite sides relative to $V_1$. The above algorithm will determine the intermediate state ${\mathbf{S}}_m={\left[\begin{array}{ll}\text{0,269901}& \text{0,469413}\end{array}\right]}^T$ and the appropriate solution structure for $W_1{W}_2=R_1{R}_2$. In other words, the intersection of integral curves gives an unsuccessful initial approximation for finding the correct ${\mathbf{S}}_m$. As a result, the algorithm leads to an inadmissible solution, shown in Fig. 2.

Fig. 2Inadmissible solution in the form of a R1R2 wave

a)

b)

Fig. 3A feasible solution in the form R1S2 wave

a)

b)

The problem lies in point 4 of the algorithm, which assumes that possible solutions that are unsuccessful from the point of view of monotonous increase in the wave velocities of the structure of the solution given by the function will arise only if there is an obligatory presence of discontinuities in the solution. The example in Fig. 1 shows that this is not so. Verify the feasibility of the solution for all possible wave configurations in the solution. In this case, a really physically correct solution, shown in Fig. 3.

Another case where the initial trial solution $W_1^{tr}{W}_2^{tr}=R_1{R}_2$ is an unsuccessful choice is shown in Fig. 4, where an enlarged fragment of the triangular diagram is shown for the same right-hand state as in Fig. 1, but with the left state ${\mathbf{S}}_l={\left(\begin{array}{ll}\text{0}& \text{0,8}\end{array}\right)}^T$. The algorithm will find a solution for $W_1^{tr}{W}_2^{tr}$, namely ${\mathbf{S}}_m={\left[\begin{array}{ll}\text{0,09734}& \text{0,793933}\end{array}\right]}^T$, and will establish that with such an intermediate state the probable structure of the solution $W_1{W}_2=R_1{S}_1{S}_2$. Indeed, this is possible because ${\mathbf{S}}_l$ and ${\mathbf{S}}_m$ are on opposite sides relatively $V_1$.

However, in the construction of the first wave $R_1{S}_1$, and then finding the intersection of the curves $S_1$ and $S_2$ in accordance with the procedures described in [3], a situation arises where the refined intermediate state ${\mathbf{S}}_m^{\left(m+1\right)}$ at the next $m+1$ step of the iteration will be on the same side relative $V_1$ to the left state. It is clear that in this case there is no sense in continuing to build a $R_1{S}_1$ wave, since the following condition must necessarily be satisfied: the left and right states must lie on opposite sides of the inflection curve $V_1$.

In fact, it is quite easy to resolve this situation. It is necessary to designate the proposed structure of the solution as incorrect and to refine it in accordance with a new approximation ${\mathbf{S}}_m$. So, for the example given in Fig. 4, after the construction of the $R_1{S}_1$-wave was not possible, after the refinement, the following structure of the solution was chosen: $W_1{W}_2=R_1{S}_2$ which turned out to be suitable, which is confirmed by Fig. 5.

Fig. 4Integral curves and Hugoniot sets for Sl=0,40,6T and Sr=0,050,9T

Fig. 5A feasible solution in the form R1S2 wave

a)

b)

3. Conclusions

The solution of the Riemann problem on an infinite domain and piecewise constant initial conditions with a single discontinuity is extremely important for practical application. A lot of laboratory experiments actually reproduce the conditions of the Riemann problem: initially the medium has homogeneous saturations, and the proportion of injected fluids remains constant during the experiment. The solution of the Riemann problem also gives information on the structure of the system of equations and can be used as a building block for obtaining solutions to problems with more complicated initial and boundary conditions. An improved algorithm for solving the Riemann problem is presented.