Shape optimization of fluid flows has a variety of applications in computational fluid mechanics, such as wing design in aerodynamics, design of bypass anastomoses in medical sciences, and pipe flows. Distributed and boundary types of Eulerian derivatives are derived from shape calculus. The distributed shape gradient in finite element discretization has better convergence and has been applied for shape design governed by Stokes, Naiver-Stokes equations, and Stokes eigenvalue problem. For approximating boundary type of Eulerian derivative, a modified boundary shape gradient formula is proposed to improve numerical accuracy based on the discrete variational outward normal derivatives. Numerical experiments are presented to verify theoretical analysis and show effectiveness of shape gradient optimization algorithms.