Convective dynamos in a rotating spherical shell feature steady zonal flows. This process is studied numerically for Prandtl numbers of 0.1 and 1, Ekman numbers in the range $E=\te{-4}$–$\te{-5}$, magnetic Prandtl number from 0.5 to 10 and Rayleigh numbers up to 100 times supercritical. The zonal flow is mainly of thermal wind origin, and minimizes the shear of the axisymmetric poloidal magnetic field lines, according to Ferraro's law of corotation. The dissipation in the interior of the fluid is mainly ohmic, while the introduction of rigid velocity boundary conditions confines viscous dissipation in the Ekman boundary layers. The root-mean-square amplitude $U_\varphi$ of the zonal flow in the spherical shell scales as $U_\varphi=(F/\Omega)^{0.5}$, F being the buoyancy flux through the shell and $\Omega$ the rotation rate. As a consequence of the corotation law, this scaling relationship is remarkably independent of the magnetic field amplitude. It does not depend on thermal, kinematic and magnetic diffusivities, owing to the large-scale and steady nature of forcing and dissipative processes. The scaling law is in agreement with the zonal-flow amplitude at the external boundary of the Earth's liquid core.