Several planetary bodies in our solar system undergo a forced libration owing to gravitational interactions with their orbital companions, leading to complex fluid motions in their metallic liquid cores or subsurface oceans. In this study, we numerically investigate flows in longitudinally librating spherical shells. We focus on the Ekman number dependencies of several shear layers when the libration frequency is less than twice of the rotation frequency and the libration amplitude is small. Time-dependent flows mainly consist of inertial waves excited at the critical latitudes due to the Ekman pumping singularities, forming conical shear layers. In particular, previous theoretical studies have proposed different scalings for the conical shear layers spawned from the critical latitudes at the inner boundary. Our numerical results favor the velocity amplitude scaling $\mathrm{O}(\varepsilon E^{1/12})$ predicted by Le Diz\`es \& Le Bars (2017) over the scaling $\mathrm{O}(\varepsilon E^{1/6})$ initially proposed by Kerswell (1995), though the Ekman numbers in our calculations are not sufficiently small to pin down this scaling. Non-linear interactions in the boundary layers drive a mean zonal flow with several geostrophic shears. Our numerical results show that geostrophic shears associated with the critical latitudes at the inner and outer boundaries exhibit the same scalings, i.e. an amplitude of $\mathrm{O}(\varepsilon^2 E^{-1/10})$ over a width of $\mathrm{O}(E^{1/5})$. Apart from the geostrophic shear associated with the critical latitude, our numerical results show that the reflection of inertial waves can induce a geostrophic shear with an amplitude of $\mathrm{O}(\varepsilon^2 E^{-1/6})$ over a width of $\mathrm{O}(E^{1/3})$.