We explore the benefits of operator splitting algorithms in the context of computational fluid dynamics. In particular, we exploit their capacity in 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 characteristics method on a fine structured grid is put forward to accurately approximate the inertial effects while continuous piecewise polynomial finite element subordinated to a coarser subdivision made of simplices is advocated for the other effects. In addition, the splitting strategy also allows to be modular and change the rheological model for the fluid in a straightforward manner.We will emphasize this flexibility by treating Newtonian flows, viscoelastic flows and multi-phase immiscible incompressible Newtonian flows based on multiple densities. 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.