<p>We propose a simple spectral element method for the 3D nonlinear spherical mean-field dynamo system, utilizing a semi-implicit discretization in time with averaging of the variable magnetic diffusivity coefficient, and a novel pointwise divergence-free approximation in space under the framework of the Helmholtz-Hodge decomposition. Based on the specifically designed basis functions, the original 3D system is decoupled into a series of independent 1D systems, which possess highly sparse matrices with no more than 18 nonzero entries in each row and multiple righthand terms, thus can be efficiently solved in parallel with the help of the fast spherical harmonics transforms. Moreover, we establish a novel spectral element approximation theory for the divergence-free vector fields in a spherical domain, and then prove rigorously the optimal error estimate of our fully discrete scheme. Finally, we present some numerical results to validate the effectiveness of our discrete scheme and main theory.</p>