Modeling of fluid flow in porous media is a pillar in geoscience applications. Previous studies have revealed that heterogeneity and fracture distribution have considerable influence on fluid flow. In this work, a numerical investigation of two-phase flow in heterogeneous fractured reservoir is presented. First, the discrete fracture model is implemented based on a hybrid-dimensional modeling approach, and an equivalent continuum approach is integrated in the model to reduce computational cost. A multilevel adaptive strategy is devised to improve the numerical robustness and efficiency. It allows up to 4-levels adaption, where the adaptive factors can be modified flexibly. Then, numerical tests are conducted to verify the the proposed method and to evaluate its performance. Different adaptive strategies with 3-levels, 4-levels and fixed time schemes are analyzed to evaluate the computational cost and convergence history. These evaluations demonstrate the merits of this method compared to the classical method. Later, the heterogeneity in permeability field, as well as initial saturation, is modeled in a layer model, where the effect of layer angle and permeability on fluid flow is investigated. A porous medium containing multiple length fractures with different distributions is simulated. The fine-scale fractures are upscaled based on the equivalent approach, while the large-scale fractures are retained. The conductivity of the rock matrix is enhanced by the upscaled fine-scale fractures. The difference of hydraulic property between homogeneous and heterogeneous situations is analyzed. It reveals that the heterogeneity may influence fluid flow and production, while these impacts are also related to fracture distribution and permeability.