The ESA mission BepiColombo will explore the planet Mercury with equipment allowing an extremely accurate tracking. While determining its orbit around Mercury, it will be possible to indirectly observe the motion of its center of mass, with an accuracy several orders of magnitude better than what is possible by radar ranging to the planet's surface. This is an opportunity to conduct a relativity experiment which will be a modern version of the traditional tests of general relativity, based upon Mercury's perihelion advance and the relativistic light propagation near the Sun. We define the mathematical methods to be used to extract from the data of the BepiColombo mission, as presently designed, the best constraints on the main post-Newtonian parameters, especially $\ensuremath{\beta},\ensuremath{\gamma}$ and the Nordtvedt parameter $\ensuremath{\eta},$ but also the dynamic oblateness of the Sun ${J}_{2\ensuremath{\bigodot}}$ and the preferred frame parameters ${\ensuremath{\alpha}}_{1},{\ensuremath{\alpha}}_{2}.$ We have performed a full cycle simulation of the BepiColombo radio science experiments, including this relativity experiment, with the purpose of assessing in a realistic (as opposed to formal) way the accuracy achievable on each parameter of interest. For $\ensuremath{\gamma}$ the best constraint can be obtained by means of a dedicated superior conjunction experiment, with a realistic accuracy $\ensuremath{\simeq}2\ifmmode\times\else\texttimes\fi{}{10}^{\ensuremath{-}6}.$ For $\ensuremath{\beta}$ the main problem is the very strong correlation with ${J}_{2\ensuremath{\bigodot}};$ if the Nordtvedt relationship $\ensuremath{\eta}=4\ensuremath{\beta}\ensuremath{-}\ensuremath{\gamma}\ensuremath{-}3$ is used, as it is legitimate in the metric theories of gravitation, a realistic accuracy of $\ensuremath{\simeq}2\ifmmode\times\else\texttimes\fi{}{10}^{\ensuremath{-}6}$ for $\ensuremath{\beta}$ and $\ensuremath{\simeq}2\ifmmode\times\else\texttimes\fi{}{10}^{\ensuremath{-}9}$ for ${J}_{2\ensuremath{\bigodot}}$ can be achieved, while $\ensuremath{\eta}$ itself is constrained within $\ensuremath{\simeq}{10}^{\ensuremath{-}5}.$ If the preferred frame parameters ${\ensuremath{\alpha}}_{1},{\ensuremath{\alpha}}_{2}$ are included in the analysis, they can be constrained within $\ensuremath{\simeq}8\ifmmode\times\else\texttimes\fi{}{10}^{\ensuremath{-}6}$ and $\ensuremath{\simeq}{10}^{\ensuremath{-}6},$ respectively, at the price of some degradation in $\ensuremath{\beta},$ ${J}_{2\ensuremath{\bigodot}}$ and $\ensuremath{\eta}.$ It is also possible to test the change with time of the gravitational constant G, but the results are severely limited because of the problems of absolute calibration of the ranging transponder, to the point that the improvement as compared with other techniques (such as lunar laser ranging) is not so important.