A numerical methodology is presented for simulating 3D multiphase flows through complex geometries on a non-body conformal Cartesian computational grid. A direct forcing implicit immersed boundary method (IBM) is used to sharply resolve complex geometries, employing the finite volume method (FVM) on a staggered grid. The fluid-fluid interface is tracked by a mass conservative sharp interface volume of fluid (VOF) method. Contact line dynamics at macroscopic length scale is simulated by imposing the apparent contact angle (static or dynamic) as a boundary condition at the three-phase contact line. The developed numerical methodology is validated for several test cases including the equilibrium shape of a droplet on flat and spherical surfaces, the temporal evolution of a droplet spreading on a flat surface. The obtained results show an excellent correspondence with those derived analytically or taken from literature. Furthermore, the present model is used to estimate, on a pore-scale, the residual oil remaining in idealized porous structures after water flooding, similar to the process used in enhanced oil recovery (EOR).