The main purpose of this work is to present a time-domain implementation of the Andrade rheology, instead of the traditional expansion in terms of a Fourier series of the tidal potential. This approach can be used in any fully three dimensional numerical simulation of the dynamics of a system of many deformable bodies. In particular, it allows large eccentricities, large mutual inclinations, and it is not limited to quasi-periodic perturbations. It can take into account an extended class of perturbations, such as chaotic motions, transient events, and resonant librations. The results are presented by means of a concrete application: the analysis of the libration of Enceladus. This is done by means of both analytic formulas in the frequency domain and direct numerical simulations. We do not a priori assume that Enceladus has a triaxial shape, the eventual triaxiality is a consequence of the satellite motion and its rheology. As a result we obtain an analytic formula for the amplitude of libration that incorporates a new correction due to the rheology. Our results provide an estimation of the amplitude of libration of the core of Enceladus as 0.6% of that of the shell. They also reproduce the observed 10 GW of tidal heat generated by Enceladus with a value of $0.17\times 10^{14}$Pa$\cdot$s for the global effective viscosity under both Maxwell and Andrade rheology.

This paper contains a review of Clairaut's theory with focus on the determination of a gravitational rigidity modulus $\gamma$ defined as $\left(\frac{C-I_o}{I_o}\right)\gamma=\frac{2}{3}\Omega^2$, where $C$ and $I_o$ are the polar and mean moment of inertia of the body and $\Omega$ is the body spin.The constant $\gamma$ is related to the static fluid Love number $k_2= \frac{3I_o G}{R^5} \frac{1}{\gamma}$, where $R$ is the body radius and $G$ is the gravitational constant. The new results are: a variational principle for $\gamma$, upper and lower bounds on the ellipticity that improve previous bounds by Chandrasekhar (1963) and a semi-empirical procedure for estimating $\gamma$ from the knowledge of $m$, $I_o$, and $R$, where $m$ is the mass of the body. The main conclusion is that for $0.2\le I_o/(mR^2)\le 0.4$ the approximation $\gamma\approx G \sqrt{ \frac{2^7}{5^5}\frac{m^5}{I_o^3}}= \gamma_I$ is a better estimate for $\gamma$ than that obtained from the Darwin-Radau equation, denoted as $\gamma_{DR}$. Moreover, within the range of applicability of the Darwin-Radau equation $0.32\le I_o/(mR^2)\le 0.4$ the relative difference between the two estimates, $|\gamma_{DR}/\gamma_I -1|$, is less than $0.05\%$.