The chaotic behavior in the rotational motion of planetary satellites is studied. A satellite is modelled as a triaxial rigid body. For a set of twelve real satellites, as well as for sets of model satellites, the full spectra of the Lyapunov characteristic exponents (LCEs) of the chaotic spatial rotation are computed numerically. The applicability of the "separatrix map approach" (Shevchenko 2000, 2002) for analytical estimation of the maximum LCEs of the rotation is studied. This approach is shown to be in a particularly good correspondence with the results of our numerical integrations in the case of a prolate axisymmetric satellite moving in an elliptic orbit. The correspondence is good in a broad range of values of the inertial parameters and the orbital eccentricity. The dependence of the LCEs on the energy (the Jacobi integral) of the system for a triaxial satellite in a circular orbit is investigated in numerical experiments. It is found that the dependence of the maximum LCEs on the energy is linear at relatively small values of the energy, with one and the same slope (but various shifts) for the majority of our sets of the values of parameters. What is more, in the case of a prolate axisymmetric satellite, the dependence seems to be universal (the same) over the studied range of the energy over a broad range of values of the inertial parameters. Upper bounds on the values of the maximum LCEs are inferred. The "energetic approach" provides a complementary method for analytical estimation of the LCEs: Evidence is given that it is useful when the energy is high.