## 3 Computational Method  
A two-dimensional model of the sphere containing phase change material was set up using ANSYS FLUENT, version 16.0, by ANSYS (Canonsburg, PA, USA). A VOF model was used for simulating the solidification and melting. The following simplifying assumptions were made for the model development:  
* Both the liquid phase and the solid phase are isotropic and homogeneous, so that constant properties can be assumed;
* Phase change material melting is considered axisymmetric;
* The two phases at the solid-liquid interface are in thermal equilibrium;
* The wall temperature of the sphere is considered constant;
* The Boussinesq approximation is adopted to simulate the convective motions inside the sphere: this consists in assuming constant density in all terms of the equations that govern natural convection, with the exception of the buoyancy term: \[(\rho-\rho_{0})g\approx-\rho_{0}\beta(T-T_{0})g\] (1) where \(\rho_{0}\) is the PCM density, \(T_{0}\) is the operating temperature and \(\beta\) is the thermal expansion coefficient. This approximation is accurate for small changes in density; in particular, the Boussinesq approximation is valid when \(\beta(T-T_{0})\ll 1\). With reference to the latter point, assuming a value of 0.001 K\({}^{-1}\) for the coefficient of thermal expansion of the A46 paraffin [32], and since the maximum temperature difference reached during the tests was 49 K, the Boussinesq approximation was appropriate for the model. Based on the previous assumptions, the model equations in spherical coordinates are [33; 34]:  
* Continuity: \[\frac{1}{r^{2}}\frac{\partial}{\partial r}\Big{(}V_{r}r^{2}\Big{)}+\frac{1}{ r\sin\theta}\frac{\partial}{\partial\theta}(V_{\theta}\sin\theta)=0\] (2)
* Polar direction momentum: \[ \begin{array}{l}\frac{\partial V_{\theta}}{\partial t}+\frac{1}{r^{2}} \frac{\partial}{\partial r}\big{(}V_{r}r^{2}V_{\theta}\big{)}+\frac{1}{r\sin \theta}\frac{\partial}{\partial\theta}(V_{\theta}\sin\theta V_{\theta})\\ =\mathrm{v}\Big{(}\frac{1}{r^{2}}\frac{\partial}{\partial r}\big{(}r^{2}\frac {\partial V_{\theta}}{\partial t}\big{)}+\frac{1}{r^{2}\sin\theta}\frac{ \partial}{\partial\theta}\big{(}\sin\frac{\partial V_{\theta}}{\partial\theta }\big{)}\Big{)}-\frac{1}{p}\frac{1}{r}\frac{\partial p}{\partial\theta}\\ -g\beta(T-T_{m})\sin\theta-\frac{V_{r}V_{0}}{r}+\
u\frac{2}{r^{2}}\frac{ \partial V_{\theta}}{\partial\theta}-\frac{v}{r^{2}\sin^{2}\theta}V_{\theta}- \frac{C(1-\alpha)^{2}}{\rho\alpha^{3}}V_{\theta}\end{array} \] (3)* Radial direction momentum: \[ \begin{array}{l}\frac{\partial V_{\mathrm{r}}}{\partial t}+\frac{1}{\mathrm{r}^{2 }}\frac{\partial}{\partial t}\big{(}V_{\mathrm{r}}\mathrm{r}^{2}V_{\mathrm{r}} \big{)}+\frac{1}{\mathrm{r}\mathrm{sin}\,\partial}\frac{\partial}{\partial \theta}\big{(}V_{\mathrm{\theta}}\mathrm{sin}\,\theta V_{\mathrm{r}}\big{)}\\ =\mathrm{v}\Big{(}\frac{1}{\mathrm{r}^{2}}\frac{\partial}{\partial r}\Big{(} \mathrm{r}^{2}\frac{\partial V_{\mathrm{r}}}{\partial t}\Big{)}+\frac{1}{ \mathrm{r}^{2}\mathrm{sin}\,\partial}\frac{\partial}{\partial\theta}\Big{(} \mathrm{sin}\,\theta\frac{\partial V_{\mathrm{r}}}{\partial\theta}\Big{)} \Big{)}-\frac{1}{\rho}\frac{\partial p}{\partial t}\\ +\mathrm{g}\mathrm{g}(\mathrm{T}-\mathrm{T}_{\mathrm{m}})\mathrm{cos}\, \theta+\frac{V_{\mathrm{r}}{}^{2}}{\mathrm{r}}-\mathrm{v}\frac{2}{\mathrm{ r}^{2}}\frac{\partial V_{\mathrm{\theta}}}{\partial\theta}-\mathrm{v}\frac{2\mathrm{ cot}\,\theta}{\mathrm{r}^{2}}V_{\mathrm{\theta}}-2\frac{\mathrm{v}}{\mathrm{r}^{2}}V_{ \mathrm{r}}\\ -\frac{\mathrm{C}(1-\alpha)^{2}}{\rho\alpha^{3}}V_{\mathrm{r}}\end{array} \] (4)
* Thermal energy equation: \[