## 3 Mathematical Formulation  
In this study, two-phase heat transfer formulation is investigated [13; 14; 15]. There are two interfaces in this study [16]. The location of the interface between air and PCM (Figure 1) is tracked by using of the volume of fluid (VOF) technique. In the VOF technique, the momentum equations along the horizontal and vertical directions with modified properties are used [17]:  
\[\begin{array}{ll}\frac{\partial(\rho u)}{\partial t}+u\frac{\partial(\rho u) }{\partial x}&+v\frac{\partial(\rho u)}{\partial y}+w\frac{\partial(\rho u)}{ \partial z}\\ &=-\frac{\partial p}{\partial x}+\rho gx+\frac{\partial}{\partial x}\Big{(} \mu\frac{\partial u}{\partial x}\Big{)}+\frac{\partial}{\partial y}\Big{(} \mu\frac{\partial u}{\partial y}\Big{)}+\frac{\partial}{\partial z}\Big{(} \mu\frac{\partial u}{\partial z}\Big{)}+\frac{10^{6}}{c^{2}+0.001}u\\ \\ \frac{\partial(\rho v)}{\partial t}+u\frac{\partial(\rho v)}{\partial x}&+v \frac{\partial(\rho v)}{\partial y}+w\frac{\partial(\rho v)}{\partial z}\\ &=-\frac{\partial p}{\partial y}+\rho gy+\frac{\partial}{\partial x}\Big{(} \mu\frac{\partial v}{\partial x}\Big{)}+\frac{\partial}{\partial y}\Big{(} \mu\frac{\partial v}{\partial y}\Big{)}+\frac{\partial}{\partial z}\Big{(} \mu\frac{\partial v}{\partial z}\Big{)}+\frac{10^{6}}{c^{2}+0.001}v\end{array} \tag{2}\]  
Due to the presence of two phases of PCM and air, modified thermophysical properties such as \(\rho,k\) and \(\mu\) can be obtained by the following equations:  
\[\rho=\beta_{air}\rho_{air}+(1-\beta_{air})\rho_{PCM} \tag{4}\]  
\[k=\beta_{air}k_{air}+(1-\beta_{air})k_{PCM} \tag{5}\]  
\[\mu=\beta_{air}\mu_{air}+(1-\beta_{air})\mu_{PCM} \tag{6}\]  
In the above equations, \(\beta_{air}\) denote the air volume fraction and obtained by the following Equation:  
\[\frac{\partial}{\partial t}[\rho_{air}\beta_{air}]+u\frac{\partial}{\partial x }[\rho_{air}\beta_{air}]+v\frac{\partial}{\partial y}[\rho_{air}\beta_{air}] +w\frac{\partial}{\partial z}[\rho_{air}\beta_{air}]=0 \tag{7}\]  
Note that the PCM volume fraction is determined with \(\beta_{PCM}=1-\beta_{air}\).  
The well-known enthalpy-porosity technique is used for the PCM phase change process [18]. In this technique, the interface between solid and liquid is not explicitly tracked. In the enthalpy-porosity approximation, the cell domain is considered as a porous medium and the amount of porosity indicates the liquid fraction \((\varepsilon)\). For \(\varepsilon=1\), the cell temperature is higher than the liquidus temperature \((T\geq T_{l})\), but the liquid fraction equal to zero indicates that the cell temperature is lower than the solidus temperature \((T\leq T_{s})\). For the cell temperature between \(T_{s}\) and \(T_{l}\), the liquid fraction \((\varepsilon)\) is between zero and one, and in this case, there is an interface between the liquid and solid phases [19]. The energy equation for problems involving phase change process is written as follow:  
\[\begin{array}{ll}\frac{\partial}{\partial t}(\rho t)+\frac{\partial}{ \partial x}(\rho u)&+\frac{\partial}{\partial y}(\rho vh)+\frac{\partial}{ \partial t}(\rho\Delta H)\\ &=\frac{\partial}{\partial x}\Big{(}k\frac{\partial T}{\partial x}\Big{)}+ \frac{\partial}{\partial y}\Big{(}k\frac{\partial T}{\partial y}\Big{)}- \frac{\partial}{\partial x}(\rho u\Delta H)-\frac{\partial}{\partial y}(\rho v \Delta H)\end{array} \tag{8}\]  
where h denote the sensible enthalpy \((h=h_{ref}+\int_{T_{ref}}^{T}c_{p}\;dT)\) and \(\Delta H\) is the latent enthalpy which is related to the amount of liquid fraction \((\Delta H=\varepsilon L)\). Note that \(L\) is the phase change enthalpy. Second-order upwind and geo reconstruct discretization methods were applied on momentum/energy and liquid fraction equations, respectively. Considering the low velocity in the phase change process, the laminar flow is considered. The convergence criteria for energy, velocity and continuity equations was selected to be \(10^{-9}\), \(10^{-4}\) and \(10^{-4}\), respectively.