The study and development of high resolution numerical approximations for the modeling and simulation of multiphase flows in petroleum reservoirs is still a challenge, from the computational viewpoint, due to the difficulties posed by some physical features such as heterogeneity and anisotropy of the medium, that are of paramount importance in this class of applications. Several methods have been proposed in the past that are based on FD (Finite-Difference), FV (Finite-Volume) or FE (Finite-Element). These methods, in their classical formulations, are of low order of approximation and suffer excessive smearing at saturation front introducing error into the numerical solution. These deficiencies can be mitigated or suppressed using highresolution methods such as the k-exact or ENO (Essentially non-Oscillatory) FV methods, which require large stencils to reconstruct high order polynomial within a control volume, resulting in an increase of the storage requirements and computational cost. On the other hand, over the last decades DG (Discontinuous Galerkin), SV (Spectral Volume), SD (Spectral Difference) and FR (Flux Reconstruction)/CPR (Correction Procedure via Reconstruction) methods were developed, which can achieve high order accuracy via a compact stencil consisting of the current cell and its immediate neighbors. In addition, the FR/CPR recovery simplified versions of nodal DG, SV and SD methods by choosing an adequate polynomial reconstruction function, whose coefficients are preprocessed and stored. The focus of this work is to investigate and to apply a very high resolution CPR method for the discretization of the saturation equation, which is generally advection-dominated and that results from the modeling of the 2-D Oil-Water displacement through porous formations. In order to suppress numerical oscillations (under/over shoots) near shocks that are typical in higher order schemes, and handing the high accuracy in smooth regions of the solution a hierarchical multi-dimensional limiting strategy (MLP) is used in the reconstruction stage. The integration in time is carried out using a third-order Runge-Kutta method. To solve the pressure equation a non-orthodox cell centered MPFA-D (Multipoint Flux Approximation-Diamond type) finite volume method is employed. In order to properly couple the MPFA-D method with the CPR formulation, it is necessary to obtain an adequate velocity reconstruction throughout the control volumes of the mesh. Because the cell-centered finite volume method naturally delivers fluxes across cell faces that belong to the primal grid, a reconstruction operator based on the lowest Raviart-Thomas interpolation functions and the Piola transformation is built, to get the complete knowledge of conservative velocity field throughout the domain. The reconstruction operator receives, as input, the density fluxes across control volume faces and returns the point-wise values of velocity anywhere within the cell. Finally, the coupling of the pressure-saturation system of equations is carried out using a classical IMPES (IMplicit Pressure Explicit Saturation) procedure. Some two-phase flow benchmark problems in one and two dimensions were analyzed and numerical and/or analytical comparisons have been used to verify the accuracy, efficiency and shock-capturing capability of the proposed methodology.FACEPEO estudo e desenvolvimento de aproximações numéricas de alta resolução para o modelagem e simulação de fluxos multifásicos em reservatórios de petróleo ainda é um desafio do ponto de vista computacional, devido às dificuldades colocadas por algumas características físicas tais como heterogeneidade e anisotropia do meio, que são de suma importância nesta classe de aplicações. Vários métodos foram propostos no passado, baseados em FD (diferenças finitas), FV (volumes finitos) ou FE (elementos finitos). Esses métodos, nas suas formulações clássicas, são de baixa ordem de aproximação e sofrem excessiva suavização na frente de saturação, introduzindo erro na solução numérica. Estas deficiências podem ser mitigadas ou suprimidas usando os métodos de alta resolução, como os métodos FV k-exact ou ENO (Essentially non- Oscillatory), que requerem estênceis grandes para reconstruir polinômios de alta ordem dentro de um volume de controle, resultando em um aumento da necessidade de armazenamento e custo computacional. Por outro lado, nas últimas décadas, métodos como DG (Discontinuous Galerkin), SV (Spectral Volume), SD (Spectral Difference) e FR (Flux Reconstruction)/CPR (Correction Procedure via Reconstruction) foram desenvolvidos, podendo alcançar alta precisão com um estêncil compacto que consiste na célula alvo e seus vizinhos imediatos. Além disso, FR/CPR recupera versões simplifixadas dos métodos nodais DG, SV e SD usando uma função de reconstrução polinomial adequada, cujos coeficientes são pré-processados e armazenados. O foco deste trabalho é investigar e aplicar o método de alta resolução CPR para a discretização da equação de saturação, que geralmente é de adveção dominante, e que resulta da modelagem do escoamento bidimensional de óleo-água em formações porosas. A fim de suprimir oscilações numéricas (under/ over shoots) perto de choques que são típicos em esquemas de alta ordem e entregar a alta precisão em regiões suaves da solução, uma estratégia limitadora multidimensional hierárquica (MLP) é usada em o estágio de reconstrução. A integração no tempo é realizada usando um método Runge-Kutta de terceira ordem. Para resolver a equação de pressão, um método dos volumes finitos não-ortodoxo o MPFA-D (Aproximação do Fluxo por Multiplos pontos-tipo Diamante) centrado na célula é empregado. Para acoplar adequadamente o método MPFA-D com a formulação do CPR é necessário obter uma reconstrução de velocidade adequada através dos volumes de controle da malha. Como o método de volumes finitos centrado na célula fornece naturalmente fluxos através das faces da célula que pertencem à malha primal, um operador de reconstrução baseado nas funções de interpolação de Raviart-Thomas de ordem mais baixa e na transformação de Piola é construído para obter o conhecimento completo do campo de velocidade conservativo em todo o domínio. O operador de reconstrução recebe, como entrada, as vazões nas faces dos volumes de controle e retorna os valores de velocidade em qualquer ponto da célula. Finalmente, o acoplamento do sistema das equações de saturação e pressão é realizado usando um procedimento clássico IMPES (IMplicit Pressure Explicit Saturation). Alguns problemas benchmark de fluxo bifásico em uma e duas dimensões são analizados e comparações numéricas e / ou analíticas foram usadas para verificar a precisão, a eficiência e a capacidade de captura de choque da metodologia proposta.El estudio y desarrollo de aproximaciones numéricas de alta resolución para el modelado y simulación de flujos multifásicos en yacimientos de petróleo sigue siendo un desafío, desde el punto de vista computacional, debido a las dificultades que presentan algunas características físicas tales como la heterogeneidad y la anisotropía del medio, que son de suma importância en esta clase de aplicaciones. Varios métodos han sido propuestos en el pasado, los cuales están basados en FD (Diferencias finitas), FV (Volumenes finitos) o FE (Elementos finitos). Estos métodos, en sus formulaciones clásicas, son de bajo orden de aproximación y sufren de excesiva dispersión en el frente de saturación introduciendo un error en la solución numérica. Esta deficiencia puede ser mitigada o suprimida usando los métodos de alta resolución como los métodos de FV k-exact o ENO (Essentially non-Oscillatory), que requieren grandes conjuntos de celdas para reconstruir polinomios de alto orden dentro de un volumen de control, lo que resulta en un aumento de la necesidad de almacenamiento y de costo computacional. Por outro lado, en las últimas décadas fueron desarrollados métodos, como, DG (Galerkin discontinuo), SV (Volumen espectral), SD (Diferencia espectral) y FR (Reconstrucción de Flujo) / CPR (Procedimiento de Corrección vía Reconstrucción), que puede lograr una precisión de alto orden a través de un conjunto de celdas compacto que consiste en la celda actual y sus vecinos inmediatos. Además, el FR / CPR recupera versiones simplificadas de los métodos nodales DG, SV y SD eligiendo una función adecuada de reconstrucción polinómica, cuyos coeficientes son preprocesados y almacenados. El objetivo de este trabajo es investigar y aplicar un método CPR de muy alta resolución para la discretización de la ecuación de saturación, que generalmente es de advección-dominante y que resulta del modelado del desplazamiento 2-D de petróleo-agua a través de formaciones porosas. Para suprimir las oscilaciones numéricas (under/ over shoots) cerca de choques que son típicos en esquemas de orden superior, y para entregar alta precisión em regiones suaves de la solución, se utiliza una estrategia de limitación multidimensional jerárquica (MLP) en la etapa de reconstrucción. La integración en el tiempo se lleva a cabo utilizando un Método de Runge-Kutta de tercer orden. Para resolver la ecuación de presión, un método de volumenes finitos no ortodoxo MPFA-D (Aproximación de flujo por multiples puntos-tipo diamante) centrado en la celda es empleado. Para poder combinar el método MPFA-D con la formulación de CPR, es necesario obtener una reconstrucción adecuada de la velocidad a través de los volúmenes de control de la malla. Devido a que el método de volumenes finitos centrado en la celda proporciona naturalmente flujos a través de las caras de las celdas que pertenecen a la malla primaria, se construye un operador de reconstrucción basado en las funciones de interpolación de Raviart-Thomas de mas bajo orden y la transformación Piola, para obtener el conocimiento completo del campo de velocidad conservativo en todo el dominio. El operador de reconstrucción recibe, como entrada, la densidad de flujos a través del volumen de control y devuelve los valores puntuales de velocidad en cualquier lugar dentro de la celda. Finalmente, el acoplamiento del sistema de ecuaciones de saturación y de presión se lleva a cabo utilizando un procedimiento clásico IMPES (IMplicit Pressure Explicit Saturation). Algunos problemas modelo para el flujo bidimensional en una y dos dimensiones se llevaron a cabo y se usaron comparaciones numéricas y / o semi-analíticas para verificar la precisión, eficiencia y capacidade de captura de discontinuidades de la metodología propuesta.