We investigate the benefits of operator splitting methods in the context of computational fluid dynamics. In particular, we exploit their capacity at handling free surface flows and a large variety of physical phenomena in a flexible way. A mathematical and computational framework is presented for the numerical simulation of free surface flows, where the operator splitting strategy allows to separate inertial effects from the other effects. The method of characteristics on a fine structured grid is put forward to accurately approximate the inertial effects while continuous piecewise polynomial finite element associated with a coarser subdivision made of simplices is advocated for the other effects. In addition, the splitting strategy also allows modularity, and in a straightforward manner rheological model change for the fluid. We will emphasize this flexibility by treating Newtonian flows, visco-elastic flows, multi-phase, and multi-density immiscible incompressible Newtonian flows. The numerical framework is thoroughly presented; the test case of the filling of a cylindrical tube with potential die swell in an extrusion process is taken as the main illustration of the advantages of operator splitting.