In this article, a new nodal discretization is proposed for two-phase Darcy flows in heterogeneous porous media. The scheme combines the Vertex Approximate Gradient (VAG) scheme for the approximation of the gradient fluxes with an Hybrid Upwind (HU) approximation of the mobility terms in the saturation equation. The discretization in space incorporates naturally nodal interface degrees of freedom (d.o.f.) allowing to capture the transmission conditions at the interface between different rock types for general capillary pressure curves. It is shown to guarantee the physical bounds for the saturation unknowns as well as a nonnegative lower bound on the capillary energy flux term. Numerical experiments on several test cases exhibit that the scheme is more robust compared with previous approaches allowing the simulation of 3D large Discrete Fracture Matrix (DFM) models.