Two new expansions of the force function of two rigid celestial bodies of finite size and arbitrary shape are obtained in Delaunay–Andoyer variables with any degree of accuracy, in the form of a partial sum of an eight dimensional Fourier series. These expansions of the force function contain products of expressions for the momenta and Stokes constants in terms of sines and cosines, whose arguments are linear combinations of the Delaunay and Andoyer angular variables. These representations of the force function are compact and convenient for applications in various problems in celestial mechanics and astrodynamics.