## 2 Material and Methods  
In this study, the physical model (shown in Figure 1) was a rectangular natural circulation loop (aspect ratio \(AR=l_{H}^{+}/l_{V}^{+}\)) containing the PCM suspension fluid. \(l_{H}^{+}\) and \(l_{V}^{+}\) were the horizontal length (width) and the vertical length (height) of the loop, respectively. The inner radius of the circular pipe in the loop was \(r_{i}^{+}\). The loop consisted of a heating and a cooling section separated by adiabatic sections. The lower horizontal section had a heating section with a length of \(l_{h}^{+}\) that provided uniformly distributed heat flux \(q_{h}^{\prime\prime}\). The upper portion of the vertical section (right) had a cooling section with length \(l_{c}^{+}\) that provided isothermal boundary \(T_{c}\). The height from the center of the cooling section to the lower horizontal section is \(l_{z}^{+}\). The geometric configuration that maximized the heat-transfer efficiency of a heating section of a rectangular thermosiphon containing PCM suspensions, obtained in our previous study [27], was used: \(AR=1\); the dimensionless length of the heating section (\(l_{h}^{*}=l_{h}^{+}/l_{H}^{+}\)) = 0.4, and the dimensionless elevation of the cooling section relative to the heating section (\(l_{z}=l_{z}^{+}/l_{c}^{+}\)) = 0.75. Consider a binary mixture system of continuous fluid (water) dispersed with small PCM particles, as shown in Figure 1 (upper-left A\(-\)A section). PCM particles can absorb/release thermal energy from/to fluid, undergoing solid-liquid phase change so that they can be solid, liquid, or both. Continuum formulation was applied to model the buoyancy-driven circulation flow of the PCM suspensions under the following assumptions. 1. PCM suspension was regarded as incompressible Newtonian fluid and in steady-state, two-dimensional, axisymmetric, and fully developed laminar flow. 2. Viscous dissipation in the loop flow was neglected. 3. All other physical properties did not vary with temperature, except that the Boussinesq approximation was used to describe working-fluid density. 4. PCM particles were evenly distributed in the loop. 5. Heat-conduction heat transfer inside the pipe wall was axisymmetric. The thermal conductivity of the pipe wall was constant. 6. There was negligible variation in the density of the PCM particles during the solid-liquid phase transition. 7. The relative velocity between PCM particles and base fluid was substantially smaller than the total velocity and could be neglected. 8. The supercooling effect of the PCMs was negligible, i.e., the melting and freezing points of the PCMs were the same. Figure 1: Schematic of rectangular thermosiphon loop and coordinate system. The following dimensionless variables were defined:  
\[r=\frac{r^{+}}{r_{i}^{+}}\,\ \text{s}= \frac{s^{+}}{l_{c}^{+}}\,\ \theta=\frac{T-T_{M}}{(q_{h}^{ \prime\prime}\ r_{i}^{+}/k_{b})}\,\ u=\frac{u^{+}}{(Ra^{\ast 1/4}\alpha_{b}/r_{i}^{+})}\]  
The assumptions and dimensionless variables given above were used to obtain the dimensionless governing equations below. Continuity equation:  
\[\frac{\partial u}{\partial s}=0. \tag{1}\]  
Momentum equation integrated over the entire loop length (\(l_{L}\)):  
\[\frac{1}{r}\frac{\partial}{\partial r}(r\frac{\partial u}{\partial r})+\frac {(Ra^{\ast 1/4}R_{L})^{3}}{\gamma_{L}l_{L}}\int_{0}^{l_{L}}\theta\phi ds=0, \tag{2}\]  
where \(\phi=1\) for \(l_{H}<s<(l_{L}/2)\); \(\phi=-1\) for \((l_{L}/2+l_{H})<s<l_{L}\); and \(\phi=0\) for \(0\leq s\leq l_{H}\) or \((l_{L}/2)<s<(l_{L}/2+l_{H})\). The correction factor for loop length \(\gamma_{L}(\geq 1)\) is embodied in the buoyancy force term of Equation (2) to specify minor losses in the loop [29].