The Laplace equation represents harmonic (i.e., both source-free and curl-free) fields. Despite the good performance of spherical harmonic series on modeling the gravitational field generated by spheroidal bodies (e.g., the Earth), the series may diverge inside the Brillouin sphere enclosing all field-generating mass. Divergence may realistically occur when determining the gravitational fields of asteroids or comets that have complex shapes, known as the Complex-boundary Value Problem (CBVP). To overcome this weakness, we propose a new numerical method based on an equivalence transformation, in which a potential-flow velocity field and a gravitational force vector field are equivalent in a mathematical sense, both referring to a harmonic vector field. The governing equation and the boundary condition of potential flow are derived here, serving as alternatives to the conventional fundamental equations, i.e., the Laplace equation and the fundamental equation of physical geodesy. Correspondingly, computational fluid dynamics (CFD) techniques are introduced as a numerical solving scheme. We apply this novel approach to the gravitational field of comet 67P/Churyumov-Gerasimenko with a singular shape. The method is validated in a closed loop simulation by comparing the result with a direct integration of Newton’s formula. It shows a good consistency between them, with a relative magnitude discrepancy at percentage level and with a maximum directional difference of 5 degrees. Moreover, the harmonicity of the potential flow’s velocity field is proved mathematically. From both theoretical and practical points of view, the new numerical method is able to overcome the divergence problem and, hence, has a good potential for solving the CBVPs.
