In this article, we consider the shape optimization problem of a body immersed in the incompressible fluid governed by Navier-Stokes equations coupling with a thermal model. Based on the continuous adjoint method, we formulate and analyze the shape optimization problem. Then we derive the structure of shape gradient for the cost functional by using the differentiability of a minimax formulation involving a Lagrange functional with the function space parametrization technique. Moreover, we present a gradient-type algorithm to the shape optimization problem. Finally, numerical examples demonstrate the feasibility and effectiveness of the proposed algorithm.