CC BY
7
Rock Flow Simulation by High-Order Quasi-Characteristics Scheme
Mikhail P. Levin <mlevin@ispras.ru> Ivannikov Institute for System Programming of RAS, 25, Alexander Solzhenitsyn st., Moscow, 109004, Russia
Abstract. A pure second-order scheme of quasi-characteristics based on a pyramidal stencil is applied to the numerical modelling of non-stationary two-phase flows through porous media with the essentially heterogeneous properties. In contrast to well-known other highresolution schemes with monotone properties, this scheme preserves a second-order approximation in regions, where discontinuities of solutions arise, as well as monotone properties of numerical solutions in those regions despite of well-known Godunov theorem. It is possible because the scheme under consideration is defined on a non-fixed stencil and is a combination of two high-order approximation scheme solutions with different dispersion properties. A special criterion according to which, one or another admissible solution is chosen, plays a key role in this scheme. A simple criterion with local character suitable for parallel computations is proposed. Some numerical results showing the efficiency of present approach in computations of two-phase flows through porous media with strongly discontinuous penetration coefficients are presented.
Keywords: Quasi-Characteristics; Two-phase Porous Media Flows; Heterogeneous Media DOI: 10.15514/ISPRAS-2018-30(5)-12
For citation: Levin M.P. Rock Flow Simulation by High-Order Quasi-Characteristics Scheme. Trudy ISP RAN/Proc. ISP RAS, vol. 30, issue 5, 2018, pp. 197-212. DOI: 10.15514/ISPRAS-2018-30(5)-12
1. Introduction
In recent fifteen years many high-resolution numerical schemes modifying Godunov scheme have been proposed (see, for instance, [1-5]). However, the problem of development of high-order schemes with monotone properties in regions near the discontinuities of solutions remains in the focus of activities for many researchers in numerical methods for partial differential equations (PDE) and in computational fluid dynamics (CFD). According to the well-known Godunov's theorem, second-order explicit monotone schemes on the fixed stencils do not exist. Up to now, two different ways to resolve this restriction are known. The first one uses the idea of lowering the approximation order in the narrow regions near the discontinuities of solution. In fact, this approach has been realized in most of the modern high resolution schemes, because they set some restrictions on the recovery functions or
limiters to provide the monotone properties of solutions in zones where the discontinuities could arise. An excellent analysis of this approach is presented in [3,5]. Therefore, most of high-resolution schemes cited above are hybrid schemes, because their approximation order is lowered in the zones near the discontinuities. Various hybrid quasi-characteristics schemes for the solution of supersonic aerodynamics problems and two-phase porous media problems were developed and considered in [6-11].
The second way consists in the construction of high-resolution schemes on the non-fixed stencils. For instance, one can apply two or more high-order schemes defined on different stencils and choose a final solution in each nodal point among admissible solutions to provide a monotone properties in regions where the discontinuities could arise. Such approach was considered in [12-14]. In these articles, various quasi-characteristics schemes of the second-order were proposed and considered. All these schemes use a combination of two second-order approximation scheme solutions having the different dispersion properties. A special criterion according to which, one or another admissible solution is chosen between two admissible solutions to provide the monotone properties near the discontinuities, plays a key role in this scheme. In [12-13], a heuristic criterion based on the third-order theoretical estimation of the average value of the governing equation operator with respect to the grid cell is proposed. Unfortunately, it has a non-local and directed character and could not be easily adopted in multidimensional case and in parallel computations. In [14], simpler local and non-directed heuristic criterions suitable for parallel computations are proposed. As is shown in [12], the quasi-characteristics schemes are more accurate than fourth-order approximation schemes in computing of the initial-value problems for the PDE of hyperbolic type, because the quasi-characteristics schemes are generalization of well-known back-ward characteristics schemes which are essentially more accurate in comparison with all other well-known numerical schemes. The reason of this consists in the naturally accurate treatment of the characteristic properties of the governing equations by the quasi-characteristics schemes in comparison with Godunov's type schemes based on the conservation laws treatment. Therefore, in recent years, various numerical schemes based on characteristics were proposed [15-19] for the solution of initial boundary value problems for reaction-diffusion equation and for correct setting of boundary conditions in decomposition of such problems.
In this article, we consider the application of a new multidimensional scheme of quasi-characteristics to the numerical simulation of two-phase flows through porous media with strongly discontinuous penetration coefficients. This scheme approximates a transport equation in the system of the porous media equations on the pyramidal stencil without any splitting. A simple criterion suitable for the selection of final solution among two admissible solutions to provide the monotone properties of the final solution without spuriuos oscillations is proposed. Numerical results for various ratio of penetration coefficients are presented. These results show
that the technique considered here could be efficiently used for the accurate numerical modelling of flows through the essentially heterogeneous porous media.
2. Governing equations, initial and boundary conditions
Let us consider the problem of a numerical simulation of two-phase flows through essentially heterogeneous porous medium with piece-wise constant absolute penetration factor. In the two-phase case, the governing equations [20] can be presented in the following form with respect to the water saturation s and the pressure p as unknown functions
,ds__д ,kkw dp, _ д ,kkw dp, _
m{dt] ~ dx{ V (1)
dx fiw Но дх dy fj,0 dy ' (2)
Here m is a porosity factor, к — k(x, y) is an absolute penetration factor of porous medium, kw = kw(s) and k0 = k0(s) are a relative penetration factors of water and oil, /iw and fj,0 are a viscosity of water and oil. Let us notice that the oil saturation s0 can be evaluated by the water saturation according to the following simple formula s0 = 1 — s.
Since the relative penetration factors kw and k0 are functions only of water saturation s, then equation (1) can be presented as follows
ds _ f к dp dk„, ч ds _ г к dp dk„, \ ds _ dt ^ mfiw dx ds ' дх ^ rn/iw dy ds ' dy
_ Ьм!-\JL (-к-_i_ d ( к Qp'il
m idx^)j,w дх> ' dy)''
Now we see that the system of governing equations (23) is of mixed type. Equation (3) is a nonlinear transport equation of hyperbolic type and the equation (2) is of elliptic type. Let us consider the transport equation (3) as a main governing equation and the elliptic equation (2) as a nonlinear restriction for coefficients of the main governing equation. Then we can apply the quasi-characteristics technique to solve the initial boundary-value problem for the transport equation and also on each time level we need to solve the boundary-value problem for elliptic equation to define the coefficients of the governing equation. In our approach, for the solution of the boundary-value problem for the elliptic equation we use the well-known five points finite difference conservative scheme and bi-conjugate gradient algorithm as in [7, 11].
Let us consider rectangle flow regions D = 0 < x < L, 0 <y < H divided into two subregions £>i = 0.2 L <x< 0.8 L, 0<y<% and D2 = D - _Db
The absolute penetration factor к in each subregion is a constant function, therefore in all region we have
(3)
k =
kDl , if (x, y) £ Di , kD2 , if (x, y) e D2 .
Initial conditions for the transport equation (3) are
six V 01 - I °-2 ' if 0<X<L, 0<y< H , S(X,y'U)~ \ 1.0 , if x = L,Q<y<H
and the boundary conditions are
H = 0 , if t> 0, y = 0,H, 0<x<L ; s(L,y,t) = 1.0 , if 0 <y<H, t > 0 .
(5)
(6)
Fig.1. Flow region
For the pressure equation (2) of the elliptic type we set a mixed Neumann and Dirichlet boundary conditions as follows
f£=0 , if 0 < x < L , y = 0,H;
p = P0,ifx = 0,0<y<H (7)
i = x = L,0<y<H.
P0 and Qo are known parameters here. The relative penetration factors of the water kw and oil ka are chosen as follows
(8)
1 f if s < 0.1 ;
kw(s) = I if 0.1 < s < 0.8
1 I 1. if s > 0.8 ;
I ' 1, if s < 0.1 ;
k0(s) = { / 0.8—5 \3 V 0.7 > ' if 0.1 < s < 0.8
, o, if s > 0.8 .
Physically, the initial boundary value problem (2-9) describes two-phase porous media flows between two horizontal wells, where the left boundary x = 0 corresponds to the production well and the right boundary x = L corresponds to the injection well.
3. Numerical scheme
In this section, we consider a non-splitting quasi-characteristics scheme on the pyramidal stencil applied to the solutions of the transport equation (3). Non-splitting scheme means that we do not use in our scheme any splitting technique for solution of the couple system of finite difference equations approximating the governing partial differential equation. It is very important in application to problems with heterogeneous coefficients, because in such problems sometimes splitting leads to the lowering of exactness of solutions. We develop this scheme with respect to the 3D transport equation written in the generalized form as follows
du , du , du , ....
m+h!Tx + b^ = h> (10)
satisfying the following initial conditions
u(0,x,y) = u0(x,y) . (11)
Here u(t, x, y) is a searching function and h = bi(t, x, y,u,^, (i = 1,2,3) and uo(x,y) are given.
In quasi-characteristics schemes [10], we approximate the governing equation written in the expanded characteristics form along some spatial grid lines I (quasi-characteristics) in (t, x. y) space as follows
(12)
Here (^г)г is a total derivative of the searching function u with respect to t along line I.
As quasi-characteristics usually are used some grid lines belonging to the considering stencil. They should lie in close vicinity with respect to the characteristics of governing equation, and sometimes can coincide with them. Now we consider a uniform, for simplicity, in each direction finite-difference grid in space (t, x, y). We denote grid steps т, hx and hy respectively. Let us consider a pyramidal stencil P1P2P3P4R in the grid space. suppose that its basement P1P2P3P4 belongs to some data layer t = t0 and vertex R belongs to the new layer t — to + т. Coordinates of the above mentioned vertices are follows: Pi(to,xo + hx,yo — hy), P2(to,xo +hx,yo +hy), Рз^о,хо — hx,yo + hy), Pi(to,xo — hx,yo — hy), R(to + т,хо,уо). Also we take into consideration m0(t0,X0, yo) a center point of the basement of the pyramid stencil and denote the nodal points corresponding to the central points of the pyramid basement ribs as
follows: m+(t0,X0 + hx,y0), m-(t0,x0 - hx,y0), n+(t0,x0,yo + hy), n-(t0,x0,y0 - hy).
We suppose that the characteristics of the transport equation going through the point R is lying inside the considering stencil and as a quasi-characteristics we can choose ribs P^R of the pyramidal stencil. Then approximating the expanded characteristic form of the governing equation along these lines, we obtain
UR-Up r/, s XR~XPjMdu\ ,
-+ K&lMfl--'
+ [(b2)PiR ~ VR /P'](ff)PiR = (h)PiR,
where i = 1,2,3,4.
(13)
Fig.2. The pyramidal stencil
According to [10] we take the following approximation of the outward derivatives at the middle layer t = to + |
du du
(^)t=to+i = {q^c + (x- X0)W + d(y - yo) , (14)
(^)t=to+i = + (y - yo)V + d(x - x0) . (15)
Here we take a center point of the middle section of our stencil (i0 + \, x0, yo) as the point C (or Ci ) and choose values W and V (or Wi and Vi ) at the middle layer io + f by the formulas
' ' Wi = W(m0) =
-2umo+um_ up2-2ti„++up3, (16)
H h^ 1 h.7 J '
v1 = v(m0) = ±[up*-2uh~++upi
v
H--TT5--1--JTS- •
We denote as W{mo) and V(rno) the finite difference operators approximating the second order derivatives of the searching function with the second order approximation error on the appropriate stencil with middle point mo.
By substitution of relations (14-17) into (13), we obtain a system of four equations with respect to ur , (f^)c, and d. Solving it we obtain ujj. In non-linear
case, we need to do three iterations on nonlinear coefficients as is usually done in the method of characteristics. We call the scheme considered above scheme I. Now we construct a second scheme of the second order approximation with different dispersive properties in comparison with scheme I. For this purpose we choose point Cu and values Wu, Vn according to the following formulas if h(x0,y0,t0) > 0 and b2(x0,yo,to) > 0 then Wu = W(P4), Vn = V(P4) , Cu — (Pix,P<iy,to + §) , if bi(x0,y0,t0) > 0 and b2(x0,yo,t0) <0 then Wu = W(P3), Vn = V(P3) ,
ClI = (P3x,P3y,to + 5) >
if bi(x0,y0,t0) < 0 and b2(xo,yo,to) <0 then Wn = W(P2), Vn = V(P2) , Cu — (P2x,p2y,to + f) , if bi(x0,y0,t0) < 0 and b2(x0,yo,to) > 0 then Wu = W(Pi), Vu = V(Pi) ,
Cu = (Plx,Ply,to + 5) ■ .
By substitution of relations (14-15) and (18) into (13), we obtain a system of four equations with respect to ur, and d. Solving it we obtain u1^. As
was mentioned above in non-linear case, it is necessary to do three iterations in evaluation of u¡(. We call this scheme scheme II.
In [12-13] for 2D case, it was shown that by choosing one of two non-monotonous admissible solutions of the second order approximation with different dispersive properties, one can construct the final solution with monotone properties. As in the papers cited before, the criterion for the choosing the final solution is based on the analysis of the average value of the governing transport equation operator evaluated on each elementary mesh cell by the high order quadrature formulas. In this criterion, the history of computations in previous grid points is taken into account and therefore is not suitable for the parallel computations and in multi-dimensional case. For 2D case, a simpler heuristic criterion based on the minimization of the rough approximation of the average value of governing operator was proposed in [14]. This criterion does not take into account the history of computations. It has a local character and it is suitable for the parallel realization.
In this paper, we construct a simple heuristic criterion as a minimal principle for the increment of searching function over the stencil in following form
(18)
UR
final _
mini=7;/7 \U1r
C0«„
m+CiMp, +С2цр2+Сзцря+С4цр4
C0+C1+C2+C3+C4
(19)
Here Ci, i = 0,1,2,3,4 are some constants to be chosen. As our numerical tests
show, the best result corresponds to the following set C0 = l, Ci = 0, ¿ = 1,2,3,4.
Thus the final solution in each grid point is chosen among two admissible solutions u\j and UjI according to the following simple minimal principle
URfinai = mini=i ii \yiR -Umo\. (20)
It is easy to see that this principle has a local character and it is very suitable for parallel computations, because it allows in principle to provide computations of searching function in each grid node independently in separate processors in computers with massive parallel processors and in computers with pipe-line processors it allows to provide the maximal loading of pipe-line.
4. Results of computations
Now let us consider some numerical results obtained by the proposed method. Computations were carried out for the following values of parameters m = 0.2 , kD2 = 1.0 ■ 10"12 m2, nw = l-10"6 N-sec- rn"2, fj,0 = 3-10"6 N -sec- to"2 , L = H = 100 m, P0= 0, Q0 = 0.69444 ■ 10"12 m2 • sec"1. Parameter kDl varies in the range from 0.50 • 10_12 m2 to 0.01 • 10"12 m2. Thus the absolute penetration in the subregion D\ is 2 to 100 times less than those in subregion D2 . Presented results correspond to the uniform grid with 61*61 nodal points in (x,y)-space.
The first series of results corresponds to fc£)l = 0.5 ■ 10_12 rn2. Fig.3 shows the isolines of water saturation s and appropriate 3D chart for time t = 400 hours. Fig.4 and 5 show the same results for t = 800 hours and t = 1200 hours respectively.
Fig.6 shows two functions characterizing the efficiency of the oil recovery process by the water drive. Line 1 corresponds to the ratio of the recovery oil to the total oil volume in initial moment t — 0 with respect to time
/[1 - s(x,y,t)]dxdy °{t) = J[l-s(x,y,0)]dxdy (21)
D
and line 2 corresponds to the function
7(i) = -j-9--(22)
describing the water content in _thg development mixture at the production well corresponding to the boundary . According to the presented results we can see that the solution of the considering problem has a wave type and the font of water wave solution is spreadjyg faster in the upper part of the flow region 2 with high penetration. Subregion 1 with low penetration plays the role of the partial obstacle and the water wave also is spre^ding^ th^e fiegi^nbut more slowly. The second series of results corresponds to Dl ' and the appropriate results
are presented on Fig.7-10. In this case, we can see that there are two shock-type water waves in the considering flow. The first wave corresponds to isolines 0.25 and 0.35 and the second corresponds to 0.45 and 0.55. According to the presented results it is easy to see that the first shock wave is spreading through the region with the low penetration, but the second wave stays near the right border of the low penetration subregion.
The third series of results corresponds to fc^ = 0.01 • 10_12 rn2 and the appropriate results are presented on fig. 11-14. In this case, the water is not spreading through the region with the low penetration and the water wave front is stopping near the right border of the low penetration subregion, which plays a role of a solid obstacle in the flow region.
Fig.3. The water saturation at t = 400 hours. ko^ = 0.50 • 10 12 m2.
Fig.4. The water saturation at t — 800 hours. ko1 = 0.50 • 10_12 TO2.
Fig.5. The water saturation at t — 1200 hours. kr>1 = 0.50 ■ 10 12 TO2.
0 500 m I500 2000 2500
t, hours
Fig.6. Characteristics of efficiency of oil recovery. Line 1 - 6{t), line 2 - 7(i). kDl = 0.50 • 10"12 TO2.
The analysis of the efficiency of the oil recovery process shows that after the moment t = 250 hours, when the water wave in the upper part of the flow region is close to the production well (boundary x = 0), the efficiency falls down and oil from the low penetration subregion and even from the high penetration subregion (0 < x < 0.2L, 0 <y< 0.5H) almost can not be developed by the water drive.
Fig.7. The water saturation at£ = 400 hours. /гд1 = 0.20 • 10 12 ш2.
Fig.8. The water saturation at £ = 800 hours. кв^ = 0.20 • 10 12 Ш2.
Fig.9. The water saturation at t = 1200 hours. ко1 = 0.20 • 10 12 Ш2.
Fig.10. Characteristics of efficiency of oil recovery. Line 1 - 0{t), line 2 - "/(t). kDl = 0.20-10"12 m2.
According to our results, we see that in the case considered in this paper, it is possible to develop only about 35 percents of oil by the usual water drive technology although 70 percents of oil is contained in the high penetration subregion. These results are in good correspondence with well-known practice.
Fig.11. The water saturation at t = 400 hours. koi = 0.01 ■ 10 12 TO2.
5. Conclusions
Our high-precision numerical quasi-characteristics technique developed for the transport equation of hyperbolic allows us to obtain solutions of complicated porous media problem with essentially heterogeneous parameters without mesh fitting procedures on rough spatial meshes. This technique can be implemented even on small computers and workstations for fast evaluation and exact modeling of oil and gas development technological processes.
Fig.13. The water saturation at t = 1200 hours. kjy1 = 0.01 ■ 10 12 m2.
Fig.14. Characteristics of efficiency of oil recovery. Line 1 - e(t), line 2 - j(t). kDl = 0.01 • 10"12 TO2.
References
[1]. B. Engquist, B. Sjogreen. High-Order Shock Capturing Methods. Computational Fluid Dynamics Review, John Willey and Sons, 1995, pp. 210-233.
[2]. E. Godlewsky, P.A. Raviart. Numerical Approximation of Hyperbolic Systems of Conservation Laws. Springer-Verlag, 1996, 524 p.
[3]. K.W. Morton. Numerical Solution of Convection-Diffusion Problems. Chapman and Hall, 1996, 384 p.
[4]. S.B. Hazra, P. Niyogi, S.K. Chakrabartty. Study in non-oscillatory schemes for shock computation using Euler equations. Computational Fluid Dynamics Journal, vol. 7, 1998, pp. 163-176.
[5]. Sh. Wo, B.M. Chen, J. Wang. A High-Order Godunov Method for One-Dimensional Convection-Diffusion-Reaction Problems. Numerical Methods for Partial Differential Equations, vol. 16, 2000, pp. 495-512.
[6]. M.P. Levin. A difference scheme of quasi-characteristics and its use to calculate supersonic gas flows. Computational Mathematics and Mathematical Physics, vol. 33, no. , 1993, pp. 113-121.
[7]. M.Yu. Zheltov, M.P. Levin. Application of the quasi-characteristics scheme for the two-phase flows through porous media. Computational Fluid Dynamics Journal, vol. 2, 1993, pp. 363-370.
[8]. M.P. Levin. Computation of 3-D supersonic flow with heat supply by explicit quasi-characteristics scheme. Computational Fluid Dynamics Journal, vol. 4, 1995, pp. 311322.
[9]. M.P. Levin, L.V. Sidorov. Hybrid modification of the scheme of the method of quasi-characteristics on a pyramidal pattern. Computational Mathematics and Mathematical Physics, vol. 35, no. 2, 1995, pp. 253- 258.
[10]. M.P. Levin. Quasi-characteristics numerical schemes. In Hyperbolic Problems: Theory, Numerics, Application, Springer, 1999, pp. 619-628.
[11]. A.I. Ibragimov, M.P. Levin, L.V. Sidorov. Numerical investigation of two-phase fluid afflux to horizontal well by quasi-characteristics scheme. Computational Fluid Dynamics Journal, vol. 8, 2000, pp. 556-560.
[12]. V.M. Borisov, Yu.V. Kurilenko, I.E. Mikhailov, E.V. Nikolaevskaya. A method of characteristics for calculation of vortex spatial supersonic stationary flows. Computing Centre of USSR Academy of Sciences, Moscow, 1988 (in Russian).
[13]. E.V. Nikolaevskaya. One class of running finite difference schemes. Computing Centre of USSR Academy of Sciences, Moscow, 1987 (in Russian).
[14]. D.Y. Kwak, M.P. Levin. High-Resolution Monotone Schemes Based on Qasi-Characteristics Technique. Numerical Methods for Partial Differential Equations, vol. 17, 2001, 262-276
[15]. S.-H. Chou, Q. Li. Characteristics-Galerkin and mixed finite element approximation of contamination by compressible nuclear waste-disposal in porous media. Numerical Methods for Partial Differential Equations, vol. 12, 1996, pp. 315-332.
[16]. H. Wang, M. Al-Lawatia, A.S. Telyakovskiy. Runge-Kutta characteristic methods for first-order linear hyperbolic equations. Numerical Methods for Partial Differential Equations, vol. 13, 1997, pp. 617-661.
[17]. H. Wang, M. Al-Lawatia, R.C. Sharpley. A characteristic domain decomposition and space-time local refinement method for first-order linear hyperbolic equations with interfaces. Numerical Methods for Partial Differential Equations, vol. 15, 1999, pp. 128.
[18]. M. Marion, A. Mollard. A multilevel characteristics method for periodic convection-dominated diffusion problems. Numerical Methods for Partial Differential Equations, vol. 16, 2000, pp. 107-132.
[19]. C.N. Dawson, M.L. Martinez-Canales. A characteristic-Galerkin approximation to a system of shallow water equations. Numerische Mathematik, vol. 86, Issue 2, 2000, pp. 239-256.
[20]. Yu.P. Zheltov. Mechanics of Oil and Gas Bearing Formation. Nedra, Moscow, 1975 (in Russian).
Численное моделирование двухфазных течений через существенно гетерогенную пористую среду схемой квазихарактеристик высокого порядка
М.П. Левин <mlevin@ispras.ru> Институт системного программирования им. В.П. Иванникова РАН, 109004, Россия, г. Москва, ул. А. Солженицына, д. 25
Аннотация. Рассматриваются вопросы численного моделирования нестационарных двухфазных потоков в пористых средах с существенно неоднородными свойствами с помощью численной схемы квазихарактеристик второго порядка аппроксимации. В отличие от известных схем высокого порядка, представленная схема имеет второй порядок аппроксимации в областях с большими градиентами решений, а также сохраняет монотонный характер решений. Это достигается за счет выбора итогового решения в каждой расчетной точке из нескольких допустимых решений с различными дисперсионными свойствами. Монотонный характер решения обеспечивается специальным критерием выбора решения, сформулированным в представленной работе. Этот критерий имеет локальный характер и удобен для параллельных вычислений. Эффективность подхода проиллюстрирована на решении задач вытеснения нефти водой в существенно неоднородных пористых пластах с коэффициентами абсолютной проницаемости, скачкообразно изменяющимися в десятки и сотни раз.
Ключевые слова: метод квазихарактеристик, параллельные вычисления, двухфазные течения, пористые гетерогенные среды.
DOI: 10.15514/ISPRAS-2018-30(5)-12
Для цитирования: Левин М.П. Численное моделирование двухфазных течений через существенно гетерогенную пористую среду схемой квазихарактеристик высокого порядка. Труды ИСП РАН, том 30, вып. 5, 2018 г., стр. 197-212 (на английском языке). DOI: 10.15514/ISPRAS-2018-30(5)-12
Список литературы
[1]. B. Engquist, B. Sjogreen. High-Order Shock Capturing Methods. Computational Fluid Dynamics Review, John Willey and Sons, 1995, pp. 210-233.
[2]. E. Godlewsky, P.A. Raviart. Numerical Approximation of Hyperbolic Systems of Conservation Laws. Springer-Verlag, 1996, 524 p.
[3]. K.W. Morton. Numerical Solution of Convection-Diffusion Problems. Chapman and Hall, 1996, 384 p.
[4]. S.B. Hazra, P. Niyogi, S.K. Chakrabartty. Study in non-oscillatory schemes for shock computation using Euler equations. Computational Fluid Dynamics Journal, vol. 7, 1998, pp. 163-176.
[5]. Sh. Wo, B.M. Chen, J. Wang. A High-Order Godunov Method for One-Dimensional Convection-Diffusion-Reaction Problems. Numerical Methods for Partial Differential Equations, vol. 16, 2000, pp. 495-512.
[6]. М. П. Левин. Разностная схема квазихарактеристик и ее применение для расчета сверхзвуковых течений газа. Ж. вычисл. матем. и матем. физ., том 33, no. 1, 1993 г., стр. 131-141.
[7]. М. Ю. Желтов, М. П. Левин, Метод квазихарактеристик в задачах механики пористых сред. Ж. вычисл. матем. и матем. физ., том 33, no.10, 1993 г., стр. 15941599.
[8]. M.P. Levin. Computation of 3-D supersonic flow with heat supply by explicit quasi-characteristics scheme. Computational Fluid Dynamics Journal, vol. 4, 1995, pp. 311322.
[9]. М. П. Левин, Л. В. Сидоров, Гибридная модификация схемы метода квазихарактеристик на пирамидальном шаблоне. Ж. вычисл. матем. и матем. физ., том 35, no.2, 1995 г., стр. 310-317.
[10]. M.P. Levin. Quasi-characteristics numerical schemes. In Hyperbolic Problems: Theory, Numerics, Application, Springer, 1999, pp. 619-628.
[11]. A.I. Ibragimov, M.P. Levin, L.V. Sidorov. Numerical investigation of two-phase fluid afflux to horizontal well by quasi-characteristics scheme. Computational Fluid Dynamics Journal, vol. 8, 2000, pp. 556-560.
[12]. В.М. Борисов, Ю.В. Куриленко, И.Е. Михайлов, Е.Л. Николаевская. Метод характеристик для расчета вихревых сверхзвуковых установившихся пространственных течений. М.: ВЦ АН СССР, 1988 г..
[13]. Е. Л. Николаевская. Об одном классе разностных схем бегущего счета. М.: ВЦ АН СССР, 1987 г..
[14]. D.Y. Kwak, M.P. Levin. High-Resolution Monotone Schemes Based on Qasi-Characteristics Technique. Numerical Methods for Partial Differential Equations, vol. 17, 2001, 262-276
[15]. S.-H. Chou, Q. Li. Characteristics-Galerkin and mixed finite element approximation of contamination by compressible nuclear waste-disposal in porous media. Numerical Methods for Partial Differential Equations, vol. 12, 1996, pp. 315-332.
[16]. H. Wang, M. Al-Lawatia, A.S. Telyakovskiy. Runge-Kutta characteristic methods for first-order linear hyperbolic equations. Numerical Methods for Partial Differential Equations, vol. 13, 1997, pp. 617-661.
[17]. H. Wang, M. Al-Lawatia, R.C. Sharpley. A characteristic domain decomposition and space-time local refinement method for first-order linear hyperbolic equations with interfaces. Numerical Methods for Partial Differential Equations, vol. 15, 1999, pp. 128.
[18]. M. Marion, A. Mollard. A multilevel characteristics method for periodic convection-dominated diffusion problems. Numerical Methods for Partial Differential Equations, vol. 16, 2000, pp. 107-132.
[19]. C.N. Dawson, M.L. Martinez-Canales. A characteristic-Galerkin approximation to a system of shallow water equations. Numerische Mathematik, vol. 86, Issue 2, 2000, pp. 239-256.
[20]. Ю.П. Желтов. Механика нефтегазоносного пласта. М.: Недра, 1975.
