A simulação numérica do fluxo de fluidos em reservatórios de petróleo naturalmente fraturados é uma ferramenta importante e com grande aplicabilidade na indústria petrolífera para o planejamento e gerenciamento da produção do reservatório. De maneira geral, a modelagem do escoamento de fluidos em meios porosos envolve a solução de um sistema de equações diferenciais parciais não-lineares que pode ser expressa por uma equação elíptica de pressão e uma equação hiperbólica de saturação, acopladas através do campo de velocidade total. Considerando que soluções semi-analíticas são possíveis apenas para casos específicos e muito simplificados, para resolver essas equações, são necessários métodos numéricos apropriados, capazes de lidar com ambientes deposicionais que apresentem grande complexidade estrutural e geométrica, o que dificulta a utilização de malhas estruturadas para representar adequadamente esses sistemas. Neste trabalho, apresentamos uma formulação numérica para simular, em domínios bidimensionais, o escoamento bifásico de óleo e água em reservatórios de petróleo naturalmente fraturados, utilizando malhas quadrilaterais não estruturadas. Para a resolução do problema difusivo da pressão, utilizou-se um Método de Volumes Finitos com Aproximação de Fluxo por Múltiplos Pontos com um estêncil de Diamante (MPFA-D), que representa uma formulação robusta e flexível, capaz de lidar com domínios fraturados altamente heterogêneos e anisotrópicos em malhas poligonais quaisquer. Para a aproximação do problema advectivo de saturação, utilizou-se o Método das Linhas de Fluxo, que é bastante eficiente devido ao desacoplamento das equações de transporte de 2-D em múltiplos problemas 1-D que são resolvidos ao longo de cada linha de corrente, reduzindo consideravelmente o custo computacional. Esses métodos foram implementados no contexto do Modelo de Fratura com Malha Híbrida (HyG), no qual, para a modelagem do escoamento bidimensional, a matriz é representada por volumes de controle de uma malha 2-D cujas arestas são alinhadas com as fraturas, que são representadas geometricamente por entidades de dimensão inferior (1-D), mas que são, posteriormente, expandidos para 2-D no espaço computacional. A formulação com o Modelo de Malha Híbrida apresentou bons resultados quando comparada com formulações similares que utilizam métodos clássicos da literatura.FACEPENumerical simulation of fluid flows in naturally fractured petroleum reservoirs is a important tool with great applicability in the oil industry for reservoir production planning and management. Modeling fluid flow can be described by a system of nonlinear partial differential equations that comprises an elliptic pressure equation and a hyperbolic saturation equation coupled through the total velocity field. Considering that semi-analytical solutions are possible only for specific and very simplified cases, to solve these equations, appropriate numerical methods are required, capable of dealing with depositional environments that present great structural and geometric complexity, which makes it difficult to use structured meshes to adequately represent these systems. In this work, we present a numerical formulation to simulate, in two-dimensional domains, the biphasic oil and water flow in naturally fractured oil reservoirs using unstructured quadrilateral meshes. For the resolution of the diffusive pressure problem, we have used a Finite-Volume Method with a Multipoint Flux Approximation with a Diamond stencil (MPFA-D), which represents a very robust and flexible formulation that is capable to deal with highly heterogeneous and anisotropic fractured domains in any polygonal grid. For the approximation of the advective saturation problem, we used the streamline-based method, which is very efficient due to the decoupling of the transport equations from 2-D into multiple 1-D problems that are solved along each streamline, considerably reducing the computational cost. These methods were implemented in the context of the Hybrid Grid Model (HyG), in which, to model the two dimensional fluid flow, the rock matrix is represented by 2- D control volumes whose edges are aligned with the fractures which are geometrically represented by entities of reduced dimension (1-D), that are later expanded to 2-D in the computational space. The formulation with the Hybrid Grid model presents very good results whenever compared with other formulations found in literature.