We present numerical results on spontaneous symmetry breaking strain localization in axisymmetric triaxial shear tests of granular materials. We simulated shear band formation using the three-dimensional distinct element method with spherical particles. We demonstrate that the local shear intensity, the angular velocity of the grains, the coordination number, and the local void ratio are correlated and any of them can be used to identify shear bands; however, the latter two are less sensitive. The calculated shear band morphologies are in good agreement with those found experimentally. We show that boundary conditions play an important role. We discuss the formation mechanism of shear bands in the light of our observations and compare the results with experiments. At large strains, with enforced symmetry, we found strain hardening.