#### Abstract

Phase change material (PCM) is particularly advantageous in developing low-carbon buildings. However, its melting efficiency is severely determined by the dynamics of the molten phase, which is greatly associated with the applied heating orientation. This work aims to describe the dynamic character in the melting process of paraffin under different heating modes. A total of four heating modes are considered for the paraffin enclosed in a square cavity, including side heating, top heating, bottom heating, and around heating. The around heating can be regarded as the combination of the first three heating modes along a single wall. A comparative study of numerical and experimental results of temperature is firstly performed for the side heating case to verify the computational model, which is then extended to study the dynamic state during the melting process of paraffin. The evolutions of the solid/liquid interface and the internal flow field show that the bottom heating can maximally activate the involvement of natural convection of the liquid paraffin and thus has the highest efficiency than the side and top heating. The contribution rate of the bottom heating is about 1.55 times over that of the top heating. This observation is useful to identify the contribution of each heating mode to the heat transfer in the paraffin and guide its practical application.

#### 1. Introduction

Stable and repeatable thermal energy storage systems, which are mainly dependent on sensible heat and latent heat storage, provide important alternatives to unrenewable fossil energy resources for balancing the supply and demand of energy and reducing the emission of CO_{2}. Compared with traditional sensible-heat-based materials, latent-heat-based phase change materials (PCMs), i.e., paraffin wax, *n*-octadecane, and tetracosane, have higher heat storage efficiency during the melting process, lower volume occupation, and relatively stable melting temperature; thus, they offer great potential in developing low-carbon buildings for efficient thermal management [1–4]. Currently, the PCMs have been integrated with the concrete wall components [5] and glazing units [6, 7] to investigate their capability in passively regulating room temperature.

However, the melting process of PCMs is extremely complex, involving heat conduction in solid and liquid phases, propagation of solid/liquid interface, thermal convective transport of liquid phase, and motion of solid phase in a fluid environment. A significant amount of experimental and numerical studies have revealed that the dynamic melting of PCM is mainly dominated by the geometric configuration of the PCM container [8] and the heating (heat source) orientations [9, 10]. For example, many of the studies have been concerned with the heat transfer of PCM in cylindrical, spherical, hexagonal, trapezoidal, rectangular, square, triangular, elliptical, annular, and hemicylindrical cavities, or even irregular nonsymmetrical containers [11–21]. Besides, the dynamic melting characteristics of enhanced PCM by nanoparticles or fin in a square cavity with partially active walls have been investigated [22, 23]. Under the specific heating temperature, the geometrical configuration of the PCM cavity can suppress or activate the natural convection of the liquid PCM, and then affect the storage or release efficiency of thermal energy. On the other hand, the relationship of the heating orientation and the gravity direction of PCM affects the melting process of PCM as well. As a widely adopted geometry, the rectangular cavity has been paid considerable attention and the melting process of PCMs in it has been investigated comprehensively when it is heated through the sidewall or the bottom wall [9, 24–26]. For the bottom heating, whose heat transfer direction is parallel but opposite to the gravity direction of PCM, the melting front of PCM becomes wavy. Comparatively, the lateral heating makes the heat transfer direction perpendicular to the gravity direction of PCM, and the melting front of PCM is an inclined curve. It has been exhibited that the thermal source from the bottom wall can maximally accelerate the melting of PCM by the induced bottom whirlpools and then shorten the melting time, compared to the side thermal source [9]. A similar effect of heating direction can be represented by changing the inclined angle of the rectangular cavity heated along one side [27]. All these studies demonstrate that the significant morphology difference of melting front is easily observed for different containers and heating directions. How to activate the natural convection effect of the liquid phase through the appropriate design of container shape and heating mode is still challenging to accelerate the phase changing efficiency and improve the heat storage capability of PCMs. Therefore, a deeper understanding of the effects of container shape and heating mode on the PCM melting mechanism is still essential to serve the design of such a highly efficient thermal management system.

In this work, we focus on studying the configuration evolution of paraffin wax (organic PCM) during the dynamic melting process in a square container under four heating modes, which is relatively rarely discussed in the literature. Paraffin wax is preferred especially for building applications, due to its wide range of melting temperature and high latent heat [28]. The PCM container is limited to a square, which is typically popular in practice to enclose the PCM for various applications. The heating mode is changed by applying heat sources in various directions to investigate how to maximize the convective behavior of liquid paraffin wax during the melting process. The paper is organized as follows. Section 2 describes the governing equations of the computational model used to simulate the heat transfer performance of paraffin wax. In Section 3, the melting test of paraffin wax under four heating modes is carried out to verify the computational model. In Section 4, the dynamic behavior of paraffin wax is investigated in detail by the verified computational model to indicate the influence of heat source direction. In Section 5, the heat transfer mechanism of paraffin for each heating mode is discussed and the heat storage efficiency is quantitatively compared, which provides a theoretical basis for the rational design of phase change heat storage system.

#### 2. Computational Model

##### 2.1. Physical Model

The calculation domain is a square cavity with side length *L* = 100 mm, as shown in Figure 1. All walls are set as nonslip conditions. In order to accurately capture the temperature gradient and velocity gradient near the nonslip hot wall, a boundary layer grid with refined cells is employed. The grid size of 1 mm is chosen to balance the computing efficiency and accuracy after the careful examination of the results of a grid refinement process, as presented in Figure 2. The whole calculation model is discretized by 15200 quadrilateral elements, including 90009 nodes. Initially, the temperature of solid paraffin is 16°C and the pressure is 1 atm. In order to improve the convergence of phase change analysis involving complex natural convection, a constant pressure point constraint is given at the left-upper point A. The heating source is a constant temperature = 70°C applied to one wall or all walls. The used PCM is paraffin wax provided by Shanghai Joule Wax Co., Ltd. The melting temperature, latent heat, and specific heat capacity were provided by the manufacturer, and the thermal conductivity was measured by the DZDR-S thermal constants analyzer in the School of Civil Engineering, Henan University of Technology. The dynamic viscosity of liquid paraffin and the coefficient of thermal expansion are chosen from the literature [26]. All material properties are listed in Table 1 for the calculation. In addition, four heating cases are considered to investigate the dynamic melting characteristic of paraffin under different boundary conditions, as shown in Table 2.

##### 2.2. Governing Equations

The melting of PCM is governed by the heat conduction equation in the solid phase and the heat conductive transport in the liquid phase described by the Navier-Stokes equations coupled with the energy equation to account for the propagated solid/liquid interface and the latent heat release during paraffin melting. In this study, the assumptions of the physical model are listed as follows: (1) The flow of fluid phase is assumed as laminar and incompressible, (2) The bulk physical properties of PCM are constant, (3) The mushy zone of PCM exists in a small temperature range, (4) The volume expansion of PCM is ignored, and (5) The buoyancy force follows the Boussinesq approximation.

The equivalent heat capacity method in the context of the fixed-domain technique is used to model phase change by the finite element software COMSOL and the involved equations include [29, 30].

Continuity equation:

Momentum equation:

Energy equation:where **u** is the velocity vector, *ρ* is the density, *t* is the time, is the Hamiltonian operator, *p* is the pressure, *μ* is the dynamic viscosity, **F** is the volume force vector, *C* is the specific heat capacity, *T* is the temperature, and *k* is the thermal conductivity.

In the phase change analysis, it is assumed that the phase change occurs in a small temperature range [*T*_{s}, *T*_{l}], and the temperature range Δ*T* = *T*_{l} − *T*_{s} generally takes 1 K, which is confirmed after careful examination of the accuracy and convergence of numerical calculation. *T*_{s} = *T*_{m} − Δ*T*/2 is the upper limit of the solid phase temperature and *T*_{l} = *T*_{m} + Δ*T*/2 is the lower limit of the liquid phase temperature. Hereafter, the subscripts *s* and *l* represent the solid and liquid phases, respectively. When *T* < *T*_{s}, it is a solid phase; when *T* > *T*_{l}, it is a liquid phase; and when *T*_{s} ≤ *T* ≤ *T*_{l} is a mixed phase (mushy zone). Then the liquid fraction *f* is defined as follows:

The volume force vector **F** in the momentum equation (2) is caused by the temperature gradient in PCM. In the solid phase, it is approximately determined by the density difference between the solid and liquid phases. In the liquid phase, it follows the popular Boussinesq approximation. In the mixed phase, it follows a linear transition.where *ρ*_{l} is the density of the liquid phase, *ρ*_{s} is the density of the solid phase, **g** is the acceleration of gravity, and *α* is the coefficient of thermal expansion.

In addition to the density variation in PCM, another problem in the phase change simulation is solved by the fixed-domain technique is the treatment of the dumping of velocity at the solid phase, where zero velocity should be ensured theoretically. In this paper, the variable viscosity method is employed. In this method, the momentum equations used keeps its original form, whereas the viscosity of the solid phase is set to be a very high value, i.e., *μ*_{s} =108*μ*_{l}. Here, *μ*_{s} and *μ*_{l} are the viscosity of the solid and liquid phases, respectively. In the mushy zone, the following rule of mixture is approximated to represent the viscosity variation:

Besides, the equivalent thermal conductivity *k* and the equivalent density *ρ* can be expressed as follows:

As the most important parameter of PCM, the equivalent specific heat *C* in energy equation (3) is used to account for phase change process, which contains the specific heat *C*_{s}, *C*_{l}, and latent heat *L*_{m}. When *T* < *T*_{s}, it is a solid phase, and the equivalent specific heat *C* is equal to the specific heat of solid phase *C*_{s}. When *T* > *T*_{l}, it is a liquid phase, and the equivalent specific heat *C* is equal to the specific heat of liquid phase *C*_{l}. When *T*_{s} ≤ *T* ≤ *T*_{l}, it is a mixed phase (mushy zone), and the specific heat is collectively determined by specific heat *C*_{s}, *C*_{l}, and latent heat *L*_{m}, as shown in the following equation:where *C*_{s} is the sensible heat of solid phase, *C*_{l} is the sensible heat of liquid phase, *L*_{m} is the latent heat of PCM, and *β* is the mass fraction, which can be expressed as follows:

##### 2.3. Numerical Approaches

In this paper, COMSOL 5.3a is used to solve the physical model discussed in the previous section. The conservation of fluxes ensured by the finite volume method makes this discretization especially advantageous when nonlinear convective terms appear as in convection dominated PCM dynamics. The convective terms are discretized using a second-order upwind scheme. The momentum and continuity equation are solved using the PIMPLE algorithm, which ensures a right pressure-velocity coupling. The under-relaxation factors for velocity, pressure, and temperature are set to 0.7, 0.3, and 0.5, respectively, these values are suitable for the particularly nonlinear problem involving natural convection. The convergence conditions for continuity, velocity, and energy equations are set to 0.0001, 0.0001, and 0.000001, respectively. A maximum temporal step of 0.1 s is imposed to ensure robustness for all the simulated conditions.

#### 3. Experimental Setup

In order to verify the numerical model, a series of melting experiments are carried out by a designed system. The experimental system includes a rectangular container of 100 × 100 × 300 mm (inner size) filled with paraffin, heat insulation plates, silicone heating plates, K-type thermocouples, a multichannel temperature acquisition device (Amber AT4108), a computer, and a digital camera, as shown in Figure 3. The rectangular container is made of 6 high-transparent plexiglass plates of the thickness of 10 mm. Besides, the container is covered by heat insulation plates to minimize energy loss during the melting process of paraffin. The heat insulation materials are white EPS foam. The thermocouples connecting the temperature acquisition device are placed in the PCM at a distance of 20 mm to collect temperate data, and the computer is to record and displace the data. The silicone heating plates were supplied by Shenzhen Hongxin Heating Technology Co., Ltd. The temperature of the heating plate was controlled by a digital temperature controller with an accuracy of ±1°C. In order to make the temperature of heating source uniform, a copper plate with thickness of 2 mm was firmly pasted to the silicone heating plate with thermal silica gel.

The melted paraffin wax is packed into the rectangular cavity and then placed at room temperature (about 16°C) at least 24 hours until the PCM is thoroughly solidified and has a uniform initial temperature. The digital temperature controller is used to control the heating temperature to maintain 70°C in the melting test, and the uncertainty of the digital temperature controller is ±1°C. Then, the experiments are carried out at room temperature.

#### 4. Results and Discussion

##### 4.1. Experimental Results

Figure 4 displays the pictures of the melting fronts of paraffin for the four heating modes when the liquid fraction of paraffin reaches 0.5. It can be found that how the configuration of the solid-liquid interface of paraffin differs when the heating condition is changed. For Case 1 in which, the square cavity is heated by constant temperature = 70°C from the right wall, the solid-liquid interface becomes inclined because the upper PCM melts faster than the bottom PCM. For Case 2, the cavity is heated by constant temperature = 70°C from the top wall. The melting front keeps almost parallel to the hot top wall. For Case 3, the cavity is heated by constant temperature = 70°C from the bottom wall. The melting front becomes irregular and exhibits a waved shape. Finally, the cavity is heated by constant temperature = 70°C around the four walls (Case 4), and it can be regarded as the combination of the former three heating modes. As a result, the melting front exhibits a more complex bell shape.

**(a)**

**(b)**

**(c)**

**(d)**

Besides, Figure 5 displays the mean local temperature evolutions overtime at six interior points 1, 2, 3, 4, 5, and 6 during the melting of paraffin for the case of heating from the right wall (Case 1). In the experiment, the heating temperature at the right wall is set as 70°C (343.15 K). Initially, the temperature at each position slowly increases up to the melting temperature of the paraffin (300 K). After that, the temperature rapidly rises until it reaches to a short plateau. Subsequently, it continues to increase and finally converges to a value very close to the heating temperature. Also, a significant difference is observed in Figure 5 for the temperature curves at different positions. At the higher position, the paraffin melts quicker and the temperature is higher. This can be attributed to the inclined solid-liquid interface, as shown in Figures 4(a) and 6(a). This means that the upper paraffin melts in advance than the lower paraffin. The thermal convective transport of the liquid phase accelerates heat transfer.

**(a)**

**(b)**

**(c)**

**(d)**

##### 4.2. Numerical Model Validation

To validate the computational model, Case 1 related to the heating from the right wall (constant temperature = 70°C) is simulated and the temperature variations at points 1–6 are recorded. It can be seen from Figure 5 that the present numerical results exhibit good agreement with the experimental data. Thus, the numerical model is adequate for simulating the concerned phase change problems with natural convection. Certainly, the small discrepancy between the numerical and experimental results is observed in Figure 5 as well, which can be explained by six possible reasons. The first one is that the desired heating temperature is difficult to be applied in the experiment. Secondly, the adiabatic condition is difficult to be kept in the experiment so the energy loss cannot be avoided thoroughly. The third one is the difference between actual and theoretical properties. The fourth one is the uncertainty in the properties measurement and temperature measurements. The fifth one is the introduced assumptions, which ignore the change of properties with respect to the temperature variable. The last one is that the three-dimensional effect is neglected in the computation. Actually, the finite distance of the front and back ends of the cavity may affect the experimental results.

##### 4.3. Numerical Results

With the verified numerical models, the detailed melting process is investigated in this section to reveal the melting mechanism for each heating mode.

###### 4.3.1. Heating from the Right Wall

Figure 7 shows the evolution of the solid-liquid interface during the melting process of paraffin when the right wall of the cavity is heated by constant temperature = 70°C. The black arrows in Figure 7 represent the velocity vector of the paraffin liquid. The heat transfers in the form of conduction at the beginning of melting, and the melting front is parallel to the hot wall at this stage. Then, the liquid paraffin starts to move upward under the drive of buoyancy force and forms a weak whirlpool, as shown in Figure 7(a). With the scouring action of horizontal flow, the solid-liquid interface at the right-top corner becomes inclined. Subsequently, the strength of the whirlpool is enhanced with time and the inclination of the melting front becomes more obvious, as shown in Figures 7(b)–7(d). During the melting process, the liquid fraction gradually increases with time. Figure 8 shows the variation of liquid fraction as a function of time. The paraffin is completely melted at *t* = 10230 s, and the melting time shortens by about 12 times over that without the consideration of natural convection (133860 s). It indicates that natural convection plays a significant role in promoting the heat storage efficiency of paraffin for the right wall heating mode. Finally, the fitted curve of the numerical results indicates that the liquid fraction can be depicted by an exponential function of time, which can be used in practice to assess the melting status of paraffin.

**(a)**

**(b)**

**(c)**

**(d)**

###### 4.3.2. Heating from the Top Wall

When the cavity is heated by constant temperature = 70°C on the top wall, the melting front of paraffin is a horizontal line and moves vertically downward. This can be indicated by Figure 9, which clearly shows the layer-distributed temperature contour and the horizontal melting front in paraffin when the liquid fraction *f* = 0.8. Such inverse temperature gradient against the gravity direction completely suppresses the generation of buoyancy force in the liquid paraffin, and thus the natural convection effect can be ignored for this heating mode. This means the heat transfers into the paraffin mainly through heat conduction. Figure 10 shows the evolution of liquid fraction in the paraffin with increasing time, and the complete melting of paraffin needs 132770 s, which is very close to the total melting time 133860 s when natural convection is discarded in the simulation. This indicates the weak natural convection effect again for the top heating. Besides, an exponent function can be used to fit the variation of liquid fraction in the paraffin in the practical application.

**(a)**

**(b)**

###### 4.3.3. Heating from the Bottom Wall

Figure 11 shows the temperature contours of the paraffin at different times when the bottom heating condition, constant temperature *T* = 70°C, is applied. It can be seen that the melting process of paraffin can be divided into four stages: (1) Heat conduction stage, (2) Stable growth stage, (3) Transition stage, and (4) Turbulent stage. Initially, the heat transfer is dominated by conduction, and no liquid move is involved. Thus, the melting front is nearly parallel to the bottom heating wall, as shown in Figure 11(a). With the thickness increase of the liquid layer, it becomes unstable and several Rayleigh-Bénard plumes are grown around the two ends of the liquid layer because of the buoyancy force caused by the temperature gradient, as shown in Figure 11(b). The number of plumes gradually increases with time, and 25 counter-rotating plumes are observed in the liquid layer at *t* = 65 s, as shown in Figure 11(c). Under the action of periodic plumes, the melting front becomes waved. This is referred as the stable growth stage. As the liquid layer further thickens, the adjacent plumes gradually merge to form a wider plume with a larger influence range, as shown in Figure 11(d)–11(e), due to the second Rayleigh-Bénard instability. In our simulation, the number of plumes has been reduced to four at *t* = 800 s after merging. As the plumes get stronger, they finally transform into a bigger turbulence and several local small turbulences [24]. The erratic motion of these turbulences makes the melting front more irregular, as shown in Figure 11(f).

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

Finally, the fitted curve of the liquid fraction in time is plotted in Figure 12 for this case. Similarly, an exponential relationship is given for practical use. The melting is completed at *t* = 4790 s, which is shortened by about 27.9 times compared to that without introducing natural convection (133860 s). This means the energy storage efficiency is greatly enhanced when the bottom heating is employed. This can be attributed to the fact that the vector direction of the temperature gradient is almost same the as the gravity direction, so the natural convection is thoroughly activated. Comparing to other single-side heating modes, the bottom heating can produce the highest melting efficiency.

###### 4.3.4. Heating from All Four Walls

The heating with constant temperature *T* = 70°C from all four walls can be regarded as the coupling of the four sing-side heating modes. The left and right wall heating modes cause the solid-liquid interface to be inclined with increasing time. The top heating makes the interface moves slowly downwards. The bottom heating makes the interface irregularly. As expected, the melting profile gradually becomes a bell shape as the melting goes on. Figure 13 shows the evolution of the solid-liquid interface with time. The total melting time is 2240 s, which is reduced by 8.9 times than that without introducing natural convection in the simulation (19910 s) for the same case.

**(a)**

**(b)**

**(c)**

**(d)**

In order to assess the contribution of each heating condition to the melting of PCM, a simple scheme is presented below. Firstly, under the four-wall heating condition, the heat flux from each heating wall is a function of the time variable and is plotted in in Figure 14. Then, the energy required for the melting of paraffin is evaluated by a simple integration as follows:where *Q*_{i} is the required energy for the *i*th heating condition, *q*_{i} is the heat flux from the heating wall, and *t*_{m} is the total melting time of paraffin for the four-wall heating condition. Here, the subscript *i* = *a*, *b*, *c,* and *d* denotes the heating from the upper wall, left wall, lower wall, and right wall, respectively.

The total melting heat storage *Q* can be expressed as follows:

Then, the contribution rate *λ*_{i} of each heating wall can be determined by a ratio of *Q*_{i} to *Q*; that is,

In this work, the contribution rates of the upper, lower, and right heating conditions are 19.2%, 29.8%, and 25.5%, respectively. As expected, the bottom heating produces the highest contribution rate, and the top heating corresponds to the smallest contribution rate. The contribution rate of the former is 1.55 times over that of the latter. This can be attributed to the strong involvement of natural convection for the bottom heating.

###### 4.3.5. Comparison of the Heating Modes

In order to analyze the effect of all heating modes on the melting efficiency of PCM, the *f*-*t* curves of cases 1–4 described above are compared in a single plot, as shown in Figure 15.

As indicated from the results in Figure 15, the total melting time of top heating, right heating, bottom heating, and all-wall heating modes are 2212.8 min, 170.5 min, 79.8 min, and 37.3 min. For the top heating mode, due to the same directions of gravity acceleration and heat source, the generation of buoyancy force, which drives the natural convective flow of molten PCM, is completely suppressed, and the longest melting time is needed to store thermal energy. As the angle between gravity acceleration and heat source increases, the effect of natural convection heat transfer is activated gradually. For the bottom heating mode, the angle between gravity acceleration and the heat source reaches its maximum, 180°. In this case, the natural convection heat transfer is completely activated, and the total melting time decreases by 96.4% compared with that of the top heating mode.

The results above indicate again that the natural convective flow driven by buoyant force plays a very important role in the melting and heat storage of phase change paraffin in a square cavity. However, due to the close correlation between the activation degree of buoyancy force and heat source direction, the thermal energy storage efficiency for different heating models exhibits a significant difference. The activation degree of natural convection heat transfer under the bottom heating condition is the highest, and this heating model provides the maximum heat storage efficiency.

#### 5. Conclusions

In this study, melting experiments and numerical simulations of paraffin enclosed in the square cavity are performed for analyzing the influence of heating modes on the transient melting characteristics of paraffin. A total of four heating modes are included: right wall heating, top heating, bottom heating, and all-wall heating. The heating temperature holds 70°C and the remaining walls keep insulated. The numerical results obtained from the equivalent heat capacity method and Boussinesq assumption are in good agreement with the experimental data. More importantly, it is found that the profile of the solid-liquid interface is dramatically influenced by the degree of involvement of natural convection, which is controlled by the heating mode. The right heating, top heating, and bottom heating cause the inclined, horizontal, and waved interfaces, respectively. For each heating condition, the liquid fraction can be described by an exponential function of time. Among all heating modes along a single wall, the bottom heating has the best thermal energy storage efficiency during the melting of the paraffin, while the top heating leads to the worst efficiency. The complete melting time for the bottom and side heating decreases by 96.4% and 92.3% than the top heating, respectively. When the square cavity is heated along all walls, the all-wall heating mode makes the flow of the liquid paraffin more complex, and the interface would become a bell shape in the melting process. The corresponding complete melting time is reduced by 98.3% than the top heating. This means that more energy can be stored if the enclosure is heated from all walls or bottom wall instead of the top wall.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The author declares that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This research was funded by the open project of the Henan Key Laboratory of Grain and Oil Storage Facility & Safety (Grant no. 2020KF-B02) and the tackled key problems in the Science and Technology Project of Henan Province, China (Grant no. 212102110203).