We study a shape optimal design problem for a forced convection flow: the steady-state Navier-Stokes equations coupled with an integro-differential thermal model. The thermal transfers are convective, diffusive and radiative with multiple reflections (model of gray bodies, radiosity equation). The inverse problem consists in minimizing a smooth cost function which depends on the solution, with respect to the domain of the equations. We prove the differentiability of the solution with respect to the domain. It follows the cost function differentiability. We introduce the adjoint state equation and obtain the exact differential of the cost function. The computational method of shape sensitivities and the optimization process are presented too.