In this paper, we propose an efficient numerical method for solving an initial boundary value problem for a coupled system of equations consisting of a nonlinear parabolic partial integro-differential equation and an elliptic equation with a nonlinear term. This problem has an important applied significance in petroleum engineering and finds application in modeling two-phase nonequilibrium fluid flows in a porous medium with a generalized nonequilibrium law. The construction of the numerical method is based on employing the finite element method in the spatial direction and the finite difference approximation to the time derivative. Newton’s method and the second-order approximation formula are applied for the treatment of nonlinear terms. The stability and convergence of the discrete scheme as well as the convergence of the iterative process is rigorously proven. Numerical tests are conducted to confirm the theoretical analysis. The constructed method is applied to study the two-phase nonequilibrium flow of an incompressible fluid in a porous medium. In addition, we present two examples of models allowing for prediction of the behavior of a fluid flow in a porous medium that are reduced to solving the nonlinear integro-differential equations studied in the paper.

The paper proposes a hybrid numerical method for solving a model problem of two-phase nonequilibrium flow of an incompressible fluid in a porous medium. This problem is relevant in the modern theory of the motion of multiphase fluids in porous media and has many applications. The studied model is based on the assumption that the relative phase permeabilities and capillary pressure depend not only on saturation, but also on its time derivative. The saturation equation in this problem refers to the type of convection-diffusion with a predominance of convection, which also includes a third-order term to account for the nonequilibrium effects. Due to the hyperbolic nature of the equation, its solution is accompanied by a number of difficulties that lead to the need for an appropriate choice of the solution method. In contrast to previous works, this paper uses a finite volume element method for solving the problem, the construction of which is based on integral balance equations, and an approximate solution is chosen from the finite element space. To discretize the problem, two different dual grids are used based on the main triangulation. In this paper, a number of a priori estimates are obtained which yields the unconditional stability of the scheme as well as its convergence with the second order. The advantages of the approach used include the local conservatism of the scheme, as well as the comparative simplicity of the software implementation of the method. These results are confirmed by a numerical test carried out on the example of a model problem.