We consider shape and topology optimization for fluids which are governed by the Navier--Stokes equations. Shapes are modelled with the help of a phase field approach and the solid body is relaxed to be a porous medium. The phase field method uses a Ginzburg--Landau functional in order to approximate a perimeter penalization. We focus on surface functionals and carefully introduce a new modelling variant, show existence of minimizers and derive first order necessary conditions. These conditions are related to classical shape derivatives by identifying the sharp interface limit with the help of formally matched asymptotic expansions. Finally, we present numerical computations based on a Cahn--Hilliard type gradient descent which demonstrate that the method can be used to solve shape optimization problems for fluids with the help of the new approach.