Heat Transfer by Natural Convection in a Square Enclosure Containing PCM Suspensions

Research on using phase change material (PCM) suspension to improve the heat transfer and energy storage capabilities of thermal systems is booming; however, there are limited studies on the application of PCM suspension in transient natural convection. In this paper, the implicit finite difference method was used to numerically investigate the transient and steady-state natural convection heat transfer in a square enclosure containing a PCM suspension. The following parameters were included in the simulation: aspect ratio of the physical model = 1, ratio of the buoyancies caused by temperature and concentration gradients = 1, Raleigh number (RaT) = 103–105, Stefan number (Ste) = 0.005–0.1, subcooling factor (Sb) = 0–1.0, and initial mass fraction (or concentration) of PCM particles (ci) = 0–0.1. The results showed that the use of a PCM suspension can effectively enhance heat transfer by natural convection. For example, when RaT = 103, Ste = 0.01, ci = 0.1, and Sb = 1, the steady-state natural convection heat transfer rate inside the square enclosure can be improved by 70% compared with that of pure water. With increasing Sb, the Nusselt number can change nonlinearly, resulting in a local optimal value.


Introduction
Heat transfer is widely used in industry and daily life, and changes in the latent heat or sensible heat of heat transfer media are used to transfer heat to the external environment. Performance improvements of a thermal system can be achieved not only by modifying the system main body itself [1], but also by either enhancing the accessories [2,3] or by changing the working fluid [4,5]. Eisapour et al. [1] numerically simulated the exergy and energy efficiencies of photovoltaic-thermal (PV-T) systems with various cooling tube configurations, including straight tubes, wavy tubes, and variable wave-amplitude wavy tubes. The influence of the geometrical parameters, the sensitivity of the operating conditions, and the working fluid effects were studied as well. Wijayanta et al. [2] experimentally investigated a double-pipe heat exchanger equipped with double-sided delta-wing (T-W) tape inserts. Liu et al. [3] numerically studied vortex flow and heat transfer enhancement in pipes with a segmented twisted element under constant wall temperature and different Reynolds numbers, ranging from 10,000 to 35,000. Muruganandam et al. [4] conducted experiments on the thermal performance enhancement of a heat exchanger by incorporating nano-particles of the natural Glay at the cold end. Ghalambaz et al. [5] aimed to improve the thermal charging performance of the conical shell-and-tube LHTES (latent heat thermal energy storage) unit by optimizing the shell shape and fin inclination angles.
Regarding the choice of working fluid, after the addition of suspended particles of the phase change material (PCM) to a working fluid, fine solid-liquid PCM particles are dispersed in a carrier or suspending working fluid (PCM slurries [6] or PCM suspensions), and the absorption and release of latent heat by PCM particles at different positions can be utilized to reduce the maximum temperature of the target facility or improve the heat transfer efficiency of the system. The manufacturing methods of PCM particles can be divided into the categories of microencapsulation technology [7] and emulsion technology [8]. In microencapsulation technology, PCM is encapsulated by polymer materials to separate the dispersed phase (PCM) and the dispersion medium (such as water and glycerin), so that the PCM particles can be suspended in the dispersion medium. Emulsion technology reduces the surface tension between the dispersed phase and the dispersion medium by adding appropriate emulsifiers, so that the PCM is uniformly suspended in the dispersion medium to form a stable emulsion.
Using a PCM suspension as a working fluid to improve heat transfer performance in forced convection in tubes, ducts, micro-channels, etc. [9,10], has been extensively displayed. Previously, studies applying the PCM suspension to natural convection were limited, and could be traced back to [11], who performed a natural convection experiment on an enclosure containing PCM particles. The heating surface was at the bottom of the enclosure, and the PCM particles were cylinders with a diameter of 3.2 mm and a height of 1.6 mm. The heat transfer coefficient obtained in this experiment was 1.9 times higher than that of pure water. Datta et al. [12] conducted a follow-up study on the experiment of [11]. The experimental model was a 150 mm × 150 mm × 150 mm cubic enclosure. The difference between these two experiments is that [12] used PCM microcapsules to replace the PCM particles used in [11]. The results showed that under low PCM concentrations (0.5% and 1%), the Nusselt number increased by approximately 80% compared with that of pure water; however, with a further increase in PCM concentration, the heat transfer efficiency gradually decreased, possibly because with increasing the PCM concentration, the uniform dispersion of the microcapsules in the fluid becomes progressively worse, and the microcapsules easily adhere to the high-or low-temperature walls, thus reducing the heat transfer efficiency. Harhira et al. [13] conducted a numerical simulation on the natural convection heat transfer from a vertical flat plate with PCM suspensions, and the results showed that the use of a working fluid containing PCM suspensions increased the Nusselt number to approximately eight times that of pure water.
Recently, PCM and PCM suspensions in natural convection have been the subject of various engineering applications, such as thermal energy storage [14,15], thermal management of buildings [16,17], heat pump and air-conditioning systems [18], electronic cooling [19][20][21], and food industry [22], and research on the use of PCM suspensions in natural convection is booming. Ho [20] investigated the heat transfer of a rectangular thermosyphon loop containing a PCM suspension using numerical simulations. The results showed that a flow regime exists where the latent heat absorption or release of the PCM suspension can effectively enhance the heat transfer, which should be defined through more extensive parametric simulations. Therefore, Ho et al. [21] investigated the effects of the loop size, the heating section proportion, and the cooling section height on the heat transfer characteristics of a rectangular thermosyphon containing a PCM suspension. The results showed that under some conditions (the aspect ratio of the rectangular loop, length of the heated section, and relative elevation of the cooled section compared to the heated section), the best heating efficiency and the cooling efficiency can be obtained.
Ghalambaz et al. [23] numerically investigated the natural convective behavior of nano encapsulated PCMs (NEPCMs) in a cavity. They showed that heat transfer enhancement is strongly related to the non-dimensional melting point. Zadeh et al. [24] numerically investigated the conjugated natural convection heat transfer and entropy generation in a square cavity composed of a porous matrix and PCM suspension with two solid blocks in it. The results illustrate that the heat transfer rates and the average Bejan number are maximized and the generated entropy is minimized when the PCM melting point is 0.5. Hajjar et al. [25] numerically studied the natural convection of a NEPCM suspension in a cavity with a hot wall with a time-periodic temperature. They indicated that the average Nusselt number in the enclosure follows a periodic variation with the same frequency as for the hot wall temperature.
Investigations on natural convection heat transfer with a PCM suspension in rectangular enclosures are gradually increasing, but most of them deal with steady-state behavior. In this context, the objective of this study is to provide insights into details of the transient buoyant flow and heat transfer in a square enclosure containing a PCM suspension, for which little information is currently available. Of particular emphasis in this study is the effects of the subcooling factor (Sb), initial concentration (c i ), Raleigh number (Ra T ), and Stefan number (Ste).

Problem Formulation
A mathematical model for the heat transfer of a square enclosure containing a PCM suspension is established, which is discretized via numerical methods and simulated using the iterative method. The physical model and boundary conditions are shown in Figure 1. The square enclosure (aspect ratio AR = H/W = 1) is filled with an aqueous fluid containing PCM particles. The upper and lower walls of the enclosure are adiabatic, the left side is a hot isothermal wall (T h ), and the right side is a cold isothermal wall (T c ). The buoyancy generated by the differences in the wall temperatures and PCM particle concentrations drive the fluid to flow in the enclosure, resulting in heat convection.
number are maximized and the generated entropy is minimized when the PCM melt point is 0.5. Hajjar et al. [25] numerically studied the natural convection of a NEP suspension in a cavity with a hot wall with a time-periodic temperature. They indica that the average Nusselt number in the enclosure follows a periodic variation with same frequency as for the hot wall temperature.
Investigations on natural convection heat transfer with a PCM suspension rectangular enclosures are gradually increasing, but most of them deal with steady-s behavior. In this context, the objective of this study is to provide insights into detail the transient buoyant flow and heat transfer in a square enclosure containing a P suspension, for which little information is currently available. Of particular emphasi this study is the effects of the subcooling factor (Sb), initial concentration ( i c ), Rale number ( T Ra ), and Stefan number (Ste).

Problem Formulation
A mathematical model for the heat transfer of a square enclosure containing a P suspension is established, which is discretized via numerical methods and simula using the iterative method. The physical model and boundary conditions are shown Figure 1. The square enclosure (aspect ratio AR = H/W = 1) is filled with an aqueous fl containing PCM particles. The upper and lower walls of the enclosure are adiabatic, left side is a hot isothermal wall (Th), and the right side is a cold isothermal wall (Tc). buoyancy generated by the differences in the wall temperatures and PCM part concentrations drive the fluid to flow in the enclosure, resulting in heat convection.  In the mixed continuum formula for the buoyancy-driven circulation flow of PCM suspension, the simplified assumptions used are as follows.
1. Suspended PCM particles are rigid spheres of a uniform size, whose size is m smaller than the characteristic dimension of flow. 2. PCM suspension is a dilute solution and is treated as a Newtonian fluid. 3. The dispersed PCM particles are neutrally buoyant in the mixture. 4. The flow of the PCM suspension is a two-dimensional laminar flow. 5. Except for the phase change process of particles, the hypothesis of local ther equilibrium between the particles and base fluid is established. 6. The melting and freezing of the PCM particles occur at a fixed temperature. 7. There are no chemical or nuclear reactions between the particles and the fluid. 8. In the initial stage, the PCM particles are uniformly suspended in the fluid, whic a two-phase fluid. 9. The interface effect of the surfactants coated on the PCM is ignored. 10. In addition to the buoyancy caused by the density differences, the thermal proper of the fluid can be regarded as constants. In the mixed continuum formula for the buoyancy-driven circulation flow of the PCM suspension, the simplified assumptions used are as follows.

1.
Suspended PCM particles are rigid spheres of a uniform size, whose size is much smaller than the characteristic dimension of flow.

2.
PCM suspension is a dilute solution and is treated as a Newtonian fluid. 3.
The dispersed PCM particles are neutrally buoyant in the mixture.

4.
The flow of the PCM suspension is a two-dimensional laminar flow.

5.
Except for the phase change process of particles, the hypothesis of local thermal equilibrium between the particles and base fluid is established. 6.
The melting and freezing of the PCM particles occur at a fixed temperature. 7.
There are no chemical or nuclear reactions between the particles and the fluid. 8.
In the initial stage, the PCM particles are uniformly suspended in the fluid, which is a two-phase fluid. 9.
The interface effect of the surfactants coated on the PCM is ignored. 10. In addition to the buoyancy caused by the density differences, the thermal properties of the fluid can be regarded as constants. 11. The density change associated with the solid-liquid phase change in the particles is negligible. 12. The thermal radiation and viscous dissipation in the flow are neglected.  13. The relative velocity between the particles and the fluid is quite low compared with the overall solution velocity, and can be ignored.
Based on the above assumptions and dimensionless variables and parameters, the governing equations can be expressed as follows: Continuity equation ∂U ∂X Energy equation Concentration equation The dimensional forms and the process of derivation of the above equations can be obtained in Ho [20]. The liquid-phase volume fraction of the PCM particles can be obtained from the following equation [26]: The effective thermal conductivity is a function of the particle concentration (c v ) and the Peclet number (Pe p ), which varies as the local shear rate is varying. Charunkyakorn et al. (1991) suggests the correlation that accounts for different values of particle Peclet number as follows: where The constants (b, m) that depend on the Peclet number are as follows: when Pe p ≤ 0.67; b = 3, m = 1.5; when 0.67 < Pe p < 250; b = 1.8, m = 0.183; when Pe p ≥ 250; b = 3, m = 0.09091. The dimensionless continuity and momentum equations can be further arranged as follows: [27] Stream function The initial conditions are set as follows: The boundary conditions are set as follows: The local and average Nusselt numbers on the hot and cold walls are defined as follows:

Numerical Method
In the grid system (as illustrated in Figure 1b), non-uniform orthogonal grid points are used for a fine resolution near the hot and cold walls. On both sides of the x-axis, a tight grid is used, and in the middle of the enclosure, a sparse grid is used; on the y-axis, the equally spaced grid is used. In this paper, the finite difference method is used to discretize the differential equations, and for the terms of time, space, and convection, the first-order backward difference, second-order central difference, and QUICK for two dimensions (QUICK-2D, [28]) are used, respectively. The vorticity on the boundaries of the walls is solved with Thom's formula. The discretization of the concentration equation on the wall uses the finite volume method. The time discretization method used is implicit discretization, and the problem-solving process is divided into two iterative loops, i.e., a time forward loop and a space iterative loop. At each time step, a spatial calculation is carried out first and must reach iterative convergence (Equation (22)), then the temporal calculation is conducted, using Equation (23) to check. If it also converges, then the iteration can advance to the next time step. The convergence criteria are defined as follows: where m represents the number of iterations and n represents the number of time steps (the iteration advances). The spatial and temporal flow field, temperature, concentration, liquid volume ratio, and vorticity are first compared with the iterative convergence criteria, which are 10 −8 for temperature and 10 −5 for other parameters. In addition, the equilibrium of the inlet and outlet heats (Nu h and Nu c ) of the high-and low-temperature walls is used to determine whether steady-state convergence is reached. The appropriate time step will depend on the parameter range. After testing, when time step ∆Fo is less than 5 × 10 −6 , the simulation results of different time steps have very little difference. Therefore, the time step used in this study is set to be 10 −6 . The spatial grid point tests are carried out for the steady states, and the average discrepancy of E Nu and E Ψ max are calculated using the changes in the average Nusselt number and the maximum value of the stream function, respectively, due to different grid-point numbers: where m1 and m2 are the different grid points.
Based on the grid test results (shown in Table 1), this study decided to use a 71 × 71 grid for the calculation.  To validate the formulation and numerical method developed in this study, the steadystate forced convection of the PCM suspension in a circular tube, as shown in Figure 2, which was reported in [29,30], is analyzed. The obtained results are qualitatively and quantitatively compared with the results from the literature. To validate the formulation and numerical method developed in this study, the steady-state forced convection of the PCM suspension in a circular tube, as shown in Figure 2, which was reported in [29,30], is analyzed. The obtained results are qualitatively and quantitatively compared with the results from the literature. The basic assumptions of this test are as follows: (1) The inlet temperature is uniform.  6) The buoyancy caused by the temperature differences is negligible. (7) The temperature field is axially symmetric.
Based on the above assumptions, the dimensionless governing equations of this problem using the methods in this paper can be expressed as follows: Energy equation: Liquid-phase volume fraction of the PCM particles: The dimensionless velocity of the fluid is completely developed: The initial and boundary conditions for this problem are: The basic assumptions of this test are as follows: (1) The inlet temperature is uniform.
(2) The flow field of the fluid is a two-dimensional, hydrodynamically fully developed laminar flow. Based on the above assumptions, the dimensionless governing equations of this problem using the methods in this paper can be expressed as follows: Energy equation: Liquid-phase volume fraction of the PCM particles: The dimensionless velocity of the fluid is completely developed: The initial and boundary conditions for this problem are: . Figure 3a shows the wall temperature distribution of the circular tube with pure water (c i = 0). The numerical simulation results of this validation model are in good agreement with the analytical solution [29] and are roughly consistent with the experimental data [29], but due to the systematic error caused by the experimental equipment, some temperatures are decreased slightly. The maximum deviation is about 20%. Figure 3b shows the influence of the Ste on the wall temperature distribution when c i = 0.1. Qualitatively, the wall temperature increases with increasing Ste. The wall temperatures obtained in this paper are slightly different from the experimental results [29], but are closer to the experimental results than the analysis results of [30]. Based on the above analysis, the predicted values of the model in this study are qualitatively consistent with the results of the literature. However, quantitatively, the differences are still significant, which may be due to the neglect of the axial heat transfer or other control variables, which is worthy of further discussion.
(31) Figure 3a shows the wall temperature distribution of the circular tube with pure water ( 0 i c = ). The numerical simulation results of this validation model are in good agreement with the analytical solution [29] and are roughly consistent with the experimental data [29], but due to the systematic error caused by the experimental equipment, some temperatures are decreased slightly. The maximum deviation is about 20%. Figure 3b shows the influence of the Ste on the wall temperature distribution when 0.1 i c = . Qualitatively, the wall temperature increases with increasing Ste. The wall temperatures obtained in this paper are slightly different from the experimental results [29], but are closer to the experimental results than the analysis results of [30]. Based on the above analysis, the predicted values of the model in this study are qualitatively consistent with the results of the literature. However, quantitatively, the differences are still significant, which may be due to the neglect of the axial heat transfer or other control variables, which is worthy of further discussion.  [29,30] shows that there was no significant difference.

Results and Discussion
The effects of Sb, i c ,

Transient Evolution
The transient evolution of the natural convection of a PCM suspension in an enclosure is affected not only by the temperature and concentration gradients, but also by the release or absorption of latent heat of the contained PCM particles. Figure 4 shows the transient evolution of the flow field, temperature distribution, and PCM solid-liquid phase distribution with time when   [29,30] shows that there was no significant difference.

Results and Discussion
The effects of Sb, c i , Ra T , and Ste on the transient and steady-state natural convection heat transfer in a square enclosure containing PCM suspension are analyzed. The governing equations used include continuity equation, momentum equations, energy equation, concentration equation, and the liquid-phase volume fraction of the PCM particles. The parameter ranges are as follows: Sb = 0-1.0, c i = 0-0.1, Ra T = 10 3 − 10 5 , and Ste = 0.005-0.1. Other parameters are set as follows: D p = 5 × 10 −4 , AR = 1, Pr = 8.5, Le = 3.3 × 10 6 , N = 1, ε = 0 (the effect of PCM supercooling is not considered), and τ = 1.

Transient Evolution
The transient evolution of the natural convection of a PCM suspension in an enclosure is affected not only by the temperature and concentration gradients, but also by the release or absorption of latent heat of the contained PCM particles. Figure 4 shows the transient evolution of the flow field, temperature distribution, and PCM solid-liquid phase distribution with time when Ra T = 10 5 , Ste = 0.05, Sb = 1, and c i = 0.05. The dotted line represents the temperature distribution, the solid line represents the flow field, the bold solid line represents the melting point, the shaded area and white area in the right figures represent the solid-phase and the liquid-phase PCM particles respectively, and the grayscale between the shaded and white areas is the phase change region. In the initial stage, the heat of the hot wall mainly transfers to the cold wall by the mechanism of heat conduction, as shown in Figure 4a. line represents the temperature distribution, the solid line represents the flow field, the bold solid line represents the melting point, the shaded area and white area in the right figures represent the solid-phase and the liquid-phase PCM particles respectively, and the grayscale between the shaded and white areas is the phase change region. In the initial stage, the heat of the hot wall mainly transfers to the cold wall by the mechanism of heat conduction, as shown in Figure 4a.  With increasing time, the distance between the melting interface and the hot wall increases, and the heat transfer mechanism changes from heat conduction to heat convection. At this moment, the velocity of PCM suspension gradually increases ( Figure  4b, max Ψ = 24.6). After the convection is enhanced, the PCM particles cross the melting interface (θ = 0;  in Figure 4c) and expand the melting range to the right (). From the figures on the right panel, the PCM inside the particles at both  is liquid. The flow turns downward as it approaches the cold wall on the right side (). During the flow process, the PCM inside the particles is gradually cooled to solid. The downward flowing With increasing time, the distance between the melting interface and the hot wall increases, and the heat transfer mechanism changes from heat conduction to heat convection. At this moment, the velocity of PCM suspension gradually increases (Figure 4b, Ψ max = 24.6). After the convection is enhanced, the PCM particles cross the melting interface (θ = 0;  With increasing time, the increases, and the heat tran convection. At this moment, th 4b, max Ψ = 24.6). After the con interface (θ = 0;  in Figure 4 figures on the right panel, the turns downward as it approa process, the PCM inside the pa PCM suspension turns to the in Figure 4c) and expand the melting range to the right (  With increasing time, the distance between the melting interface and the hot wall increases, and the heat transfer mechanism changes from heat conduction to heat convection. At this moment, the velocity of PCM suspension gradually increases ( Figure  4b  With increasing time, the distance between the melting interface increases, and the heat transfer mechanism changes from heat co convection. At this moment, the velocity of PCM suspension gradually 4b, max Ψ = 24.6). After the convection is enhanced, the PCM particles interface (θ = 0;  in Figure 4c) and expand the melting range to the ri figures on the right panel, the PCM inside the particles at both  is turns downward as it approaches the cold wall on the right side (). process, the PCM inside the particles is gradually cooled to solid. The do PCM suspension turns to the left under the influence of the wall belo With increasing time, the distance between the melting interface and the hot wall increases, and the heat transfer mechanism changes from heat conduction to heat convection. At this moment, the velocity of PCM suspension gradually increases ( Figure  4b, max Ψ = 24.6). After the convection is enhanced, the PCM particles cross the melting interface (θ = 0;  in Figure 4c) and expand the melting range to the right (). From the figures on the right panel, the PCM inside the particles at both  is liquid. The flow turns downward as it approaches the cold wall on the right side (). During the flow process, the PCM inside the particles is gradually cooled to solid. The downward flowing PCM suspension turns to the left under the influence of the wall below (), then it is is liquid. The flow turns downward as it approaches the cold wall on the right side ( With increasing time, the distance between the melting interface and the hot wall increases, and the heat transfer mechanism changes from heat conduction to heat convection. At this moment, the velocity of PCM suspension gradually increases ( Figure  4b, max Ψ = 24.6). After the convection is enhanced, the PCM particles cross the melting interface (θ = 0;  in Figure 4c) and expand the melting range to the right (). From the figures on the right panel, the PCM inside the particles at both  is liquid. The flow turns downward as it approaches the cold wall on the right side (). During the flow process, the PCM inside the particles is gradually cooled to solid. The downward flowing PCM suspension turns to the left under the influence of the wall below (), then it is ). During the flow process, the PCM inside the particles is gradually cooled to solid. The downward flowing PCM suspension turns to the left under the influence of the wall below ( With increasing time, the distance between the melting interface and the hot wall increases, and the heat transfer mechanism changes from heat conduction to heat convection. At this moment, the velocity of PCM suspension gradually increases ( Figure  4b, max Ψ = 24.6). After the convection is enhanced, the PCM particles cross the melting interface (θ = 0;  in Figure 4c) and expand the melting range to the right (). From the figures on the right panel, the PCM inside the particles at both  is liquid. The flow turns downward as it approaches the cold wall on the right side (). During the flow process, the PCM inside the particles is gradually cooled to solid. The downward flowing PCM suspension turns to the left under the influence of the wall below (), then it is ), then it is directed upward by the thermal buoyancy near the hot wall on the left (   11 of 18 directed upward by the thermal buoyancy near the hot wall on the left (). Finally, it returns to the top region () to form the main circulation.
When Fo = 0.05 (Figure 4c), there is an inner circulation inside the above-mentioned main circulation (symbol (A)), and the right figure of Figure 4c shows that this small circulation is formed by solid PCM particles. As time passes, when Fo = 0.1 (Figure 4d), this internal small circulation becomes two small circulations (symbols (B) and (C)), and the right figure of Figure 4d shows that the left circulation (B) is formed by solid PCM particles, while the right circulation (C) is formed by liquid PCM particles, and there is a solid PCM region inside the small circulation (C). When Fo = 0.2, as shown in the right figure of Figure 4e, the aforementioned solid PCM region is melted, and a small circulation (D) is fully formed by liquid PCM particles. When the circulation crosses the melting interface, the PCM particles in the working fluid could be in different phases before and after the crossing. Figure 4f (steady-state) shows that the internal circulation at (E) is full of liquid PCM particles. When the fluid flows to the right and contacts the lowtemperature region, the particles dissipate heat. When the fluid flows back to the left region at (F), the PCM particles have solidified.

Steady-State Results
In this section, analyses are performed on the effects of Sb, i c ,

Steady-State Results
In this section, analyses are perform steady-state results. It can be seen from increases; however, with increasing Sb, N local optimum Nu . When When Fo = 0.05 (Figure 4c), there is an inner circulation inside the above-mentioned main circulation (symbol (A)), and the right figure of Figure 4c shows that this small circulation is formed by solid PCM particles. As time passes, when Fo = 0.1 (Figure 4d), this internal small circulation becomes two small circulations (symbols (B) and (C)), and the right figure of Figure 4d shows that the left circulation (B) is formed by solid PCM particles, while the right circulation (C) is formed by liquid PCM particles, and there is a solid PCM region inside the small circulation (C). When Fo = 0.2, as shown in the right figure of Figure 4e, the aforementioned solid PCM region is melted, and a small circulation (D) is fully formed by liquid PCM particles. When the circulation crosses the melting interface, the PCM particles in the working fluid could be in different phases before and after the crossing. Figure 4f (steady-state) shows that the internal circulation at (E) is full of liquid PCM particles. When the fluid flows to the right and contacts the low-temperature region, the particles dissipate heat. When the fluid flows back to the left region at (F), the PCM particles have solidified. Figure 5 illustrates the parameter effects (Sb, c i , Ra T , and Ste) on the average Nusselt numbers of hot and cold walls. Figure Figure 5c shows the changes in Nu under different Ra T s. The increases in Ra T can result in a strong convection effect, which further increases the heat transfer rate. As shown in the figure, the time required to reach equilibrium under a large Ra T is shorter than that in other cases, and the higher the Ra T , the higher the heat transfer rate is. When Ra T is increased from 10 3 to 10 5 , the steady-state Nu is increased from 1.05 to 5.0, showing substantial growth. However, cases with high Ra T are more susceptible to the influence of other parameters regarding the generation of oscillatory natural convection. It can be seen from Figure 5d that with decreasing Ste, Nu increases, and when Ste = 0.01, the steady-state Nu is 1.9. According to above analysis, under the condition of high c i , high Ra T , and low Ste, the transient and the steady-state values of Nu are high, and the time required to reach steady-state is shortened.

Steady-State Results
In this section, analyses are performed on the effects of Sb, c i , Ra T , and Ste on the steady-state results. It can be seen from Figure 6 that with increasing c i or Ra T , Nu h increases; however, with increasing Sb, Nu h increases first then decreases, and there is a local optimum Nu. When Ra T = 10 3 (Figure 6a), regardless of c i , the optimal heat transfer (Nu h ) always occurs at Sb = 0.8; when Ra T = 10 4 (Figure 6b, the optimum heat transfer characteristics occur at Sb = 0.5; and when Ra T = 10 5 (Figure 6c), the optimal heat transfer characteristics occur at Sb = 0.5. In summary, within the discussed parameter ranges, when Ra T = 10 3 , the maximum Nu h is 1.20, which occurs when c i = 0.1 and Sb = 0.8; when Ra T = 10 4 , the maximum Nu h is 2.89, which occurs when c i = 0.1 and Sb = 0.5; and when Ra T = 10 5 , the maximum Nu h is 5.13, which occurs when c i = 0.1 and Sb = 0.5.     Nusselt number also decreases, which is consistent with the observations by Katz [11]. An increase in T Ra means that the convection effect of the PCM particles in the system is increased. As a result of the enhanced convection, the time for the PCM particles to  Table 2 shows the effects of Ra T and Ste. According to the definition of Ste, a low Ste means that more latent heat than sensible heat is contained in the fluid, and with increasing latent heat, the overall heat transfer characteristics can be improved. According to the data in the table, the lower the Ste, the higher the heat transfer characteristics. When Ra T = 10 3 and Ste = 0.01, the maximum increase in the Nusselt number of the high-temperature wall can be approximately 70% higher than that of pure water (Nu h /Nu w = 1.699, the highest value within the discussed parameter ranges). As Ra T increases, the Nusselt number also decreases, which is consistent with the observations by Katz [11]. An increase in Ra T means that the convection effect of the PCM particles in the system is increased. As a result of the enhanced convection, the time for the PCM particles to experience a phase change in the phase change region decreases, which reduces the release and absorption of latent heat, thereby reducing the overall heat transfer capability. As Ra T increases, Ste increases, which can induce oscillatory natural convection; in other words, low latent heat can induce oscillatory natural convection under high Ra T . This phenomenon is not the focus of this study, but it is worthy of subsequent in-depth investigation.

Conclusions
Research on natural convection heat transfer with PCM suspension in rectangular enclosures are gradually increasing; however, there are limited studies on the transient behavior. In this context, the objective of this study focused on the transient and steady-state natural convection of PCM suspensions in a square enclosure, mainly influenced by the Raleigh number (Ra T ), Stefan number (Ste), subcooling factor (Sb), and initial concentration (c i ). The following parameters were included: aspect ratio of the physical model = 1, ratio of the buoyancies caused by temperature and concentration gradients = 1, Ra T = 10 3 -10 5 , Ste = 0.005-0.1, Sb = 0-1.0, and c i = 0-0.1. The results can be summarized as follows: 1.
Under the condition of high c i , high Ra T , and low Ste, the transient and the steadystate values of the Nusselt number (Nu) are high, and the time required to reach steady-state is decreased.

2.
With increasing Sb, Nu increases first then decreases, and there is a local optimum N u . Within the discussed parameter ranges, when Ra T = 10 3 , the maximum Nu h is 1.20, which occurs when c i = 0.1 and Sb = 0.8; when Ra T = 10 4 , the maximum Nu h is 2.89, which occurs when c i = 0.1 and Sb = 0.5; and when Ra T = 10 5 , the maximum Nu h is 5.13, which occurs when c i = 0.1 and Sb = 0.5.

3.
As Ra T increases, Nu h decreases; moreover, as Ra T increases, oscillatory natural convection is induced by higher Stes; in other words, the low latent heat of the PCM can induce oscillatory natural convection under high Ra T . This phenomenon of oscillatory natural convection has not yet been mentioned in the literature.
The modeling and simulation results in this study are limited to the PCM used, the geometric configuration, and the chosen parameter ranges. An in-depth investigation on the behavior of PCM slurries or suspensions in practical applications is noticeable. This is not what we are exploring in this study, but it is worthy of further consideration.