Reynolds-Averaged Navier-Stokes (RANS) equations are widely used in engineering for turbulent flow simulations due to the significantly reduced computational effort to discretize them compared classic Navier-Stokes equations. The unknows of the RANS equations are the averaged velocity, the pressure fields respectively and the Reynolds stress tensor (RST). The latter is a symmetric tensor that needs to be modeled to close the RANS equations. Classically, this is performed by deducing additional partial differential equations obtained heuristically through physical arguments. While these models have decades of history and they are well-established, they can be highly inaccurate for some types of flow. Recently, several studies have used machine learning techniques to learn the mapping between the averaged fields and the RST to increase the accuracy of RANS models [1]. In this talk, we present a data-driven model based on neural networks that predicts the divergence of the RST and guarantees both Galilean and frame-reference rotation invariances by construction. Numerical experiments are presented to show the improvements of our approach compared to standard RANS models.