A Finite-Difference Scheme for Three-Dimensional Incompressible Flows in Spherical Coordinates

In geophysics and astronomy, solving the Navier-Stokes equations in spherical coordinates is a desirable option owing to the simmetry properties of boundaries and discretization. Unfortunately, these equations have two singularities, at the centre and along the polar axis, that make their integration exceedingly difficult. In this paper, a second-order finite-difference scheme for 3D, incompressible flows in spherical coordinates is presented. Thanks to a staggered mesh and a change of variables, the singular boundaries become trivial boundary conditions for the modified variables. The integration is performed by a fractional-step method with implicit viscous- and explicit nonlinear-terms. The elliptic equation for incompressibility benefits from a direct solver thus obtaining a free divergence velocity field within round-off error. The algorithm is efficient and flexible allowing any mesh distribution in two of the three spatial directions. The method has been validated by a Hill vortex crossing the singular boundaries and comparing the results with simulations in Cartesian coordinates and the theory. The scheme has been used also for Rayleigh-Bénard convection within spherical shells and the results compared with the literature.