2.1. Description of physical problem
Modeling creates a simplified representation of a problem that can generally be defined in mathematical model. This model will simulate the physical phenomena involved and transcribed in mathematical or physical form to predict the particular constraints of a given situation. Regarding latent thermal storage, many numerical simulations do not take into account some physical phenomena, including the often ignored PCMs hysteresis.
Indeed, the modeling study, taking this phenomenon into account, makes the thermal performance evaluation of PCMs more objective and reliable in simulation studies. It also allows a better understanding of the phenomena surrounding the exploits of these materials, and it allows more efficiency during the application. Therefore, the main objective of this section is to provide a simulation method with MATLAB program for the different physical phenomena involved in the thermal modeling of the building to adapt its energy behavior, taking into account the environmental conditions, through the following steps:
- Discussion on the numerical method used to model the thermal behavior of the PCM,
- Study of the temperature evolution through a complex wall containing an PCM layer by coupling the conduction, convection and radiation models; where the hysteresis is taken into account in the modeling of PCM,
- Model validation using the dynamic simulation program EnergyPlus Software.
Figure 2 shows the physical problem treated in the modeling process, the structure considered in the study is a building wall consisting of a layer of concrete with a thickness of 0.15 m, covered with a thin layer of 0.02 m thick PCM on the inner face. Table 1 reviews the thermophysical properties of the materials used in the simulation. It should be noted that the indoor air temperature was considered constant equal to 24°C. The same applies to the internal and external convection coefficients, estimated respectively at 3.079 and 11 W/m2 K. As for other factors related to weather conditions such as solar radiation, sky temperature and clarity, etc., to establish the boundary conditions on the outer face, they are implicitly calculated in the algorithm according to the formulas in the literature, which was mentioned in detail in [11]. This makes the study results very useful for realistic practical applications of PCMs in building structures.
Table 1
Thermophysical properties of concrete and PCM.
Material | Concrete | PCM |
---|
Thermal conductivity (W/m K) | 0.7330 | 0.7264 |
Density (kg/m3) | 2315 | 1601.85 |
Specific heat capacity (J/kg K) | 800 | 836.86 |
Latent heat of fusion (kJ/kg) | - | 137.40 |
Melting temperature (°C) | - | 22 (or 20 in case of solidification) |
Thickness (m) | 0.15 | 0.02 |
2.2. Mathematical formulation
In this study, the enthalpy method was chosen for mathematical formulation, which in turn is one of the most popular methods for modeling nonlinear phase change behavior in PCMs, was first introduced by Eyres et al. [12] to deal with changes in thermal properties with respect to temperature. In the enthalpy method, the total amount of energy required during the phase-change process is dealt with by including both latent and specific heat in the enthalpy term in the governing equation (Eq. (1)).
\(\rho \frac{{\partial h}}{{\partial t}}=\frac{\partial }{{\partial x}}\left( {k\frac{{\partial T}}{{\partial x}}} \right)............................................................Eq.(1)\)
Regarding the assumptions established in this study, a one-dimensional heat transfer model is assumed with a totally implicit time step. The spatial discretization is based on the finite volume method and the harmonic mean suggested by Patankar [13] for the conductivity of materials is used. Figure 3 illustrates typical grid points of a wall layer system using the finite volume method.
Convergence is declared using the following relation:\(Erreur=\left| {\frac{{\sum {\left( {T - {T_{new}}} \right)} }}{{\sum {\left( {{T_{new}}} \right)} }}} \right|\)
Where T is the temperature vector of the nodes of the previous iteration and Tnew is the new results.
Using Fig. 3(b), a fully implicit approximation of the control volume is used. Therefore, for internal nodes, Eq. (1) can be discretized as follows:
\(\rho \frac{{h_{p}^{{n+1}} - h_{p}^{n}}}{{\Delta t}}={k_w}\frac{{T_{w}^{{n+1}} - T_{p}^{{n+1}}}}{{\Delta X \times \delta {X_w}}}+{k_e}\frac{{T_{e}^{{n+1}} - T_{p}^{{n+1}}}}{{\Delta X \times \delta {X_e}}}..................................................Eq.(2)\)
With:
n: indicate a time step, w: west node, p: point node, e: east node.
By putting the terms together and rearranging them, Eq. (2) becomes:
\(h_{p}^{{n+1}}=h_{p}^{n}+{a_w} \times T_{w}^{{n+1}}+{a_p} \times T_{p}^{{n+1}}+{a_e} \times T_{e}^{{n+1}}.....................................Eq.(3)\)
With:
\({a_w}=\frac{{{k_w} \times \Delta t}}{{\rho \times \Delta X \times \delta {X_w}}},\) \({a_e}=\frac{{{k_e} \times \Delta t}}{{\rho \times \Delta X \times \delta {X_e}}},\) \({a_p}= - \left( {{a_w}+{a_e}} \right)\)
Using the method proposed by Patankar [13] to linearize Eq. (3), because hpn+1and T n+1 are unknown at this time step. The term hpn+1 can be developed using the first order approximation of the Taylor series as follows:
\(h_{p}^{{n+1}}=h_{p}^{n}+C{(T)^{n+1,m}} \times \left( {T_{p}^{{n+1}} - T_{p}^{n}} \right)...................................................Eq.(4)\)
At the start of the simulation, the temperature fields are based on approximate values. Therefore, Cp (provisional values that will be updated at the start of the iteration) can be found using:
\(C\left( T \right)=\left\{ \begin{gathered} {C_S}..................................................T \leqslant {T_m} - \varepsilon \hfill \\ \frac{{{C_S}+{C_l}}}{2}+\frac{L}{{2 \times \varepsilon }}..............................{T_m} - \varepsilon <T<{T_m}+\varepsilon .......................Eq.(5) \hfill \\ {C_l}...................................................T \geqslant {T_m}+\varepsilon \hfill \\ \end{gathered} \right.\)
By substituting the function hpn+1 in Eq. (3) and after simplifying and rearranging the terms, we obtain the following estimated linear equation:
\(a_{w}^{{n+1,m}} \times T_{w}^{{n+1,m+1}}+a_{p}^{{n+1,m}} \times T_{p}^{{n+1,m+1}}+a_{e}^{{n+1,m}} \times T_{e}^{{n+1,m+1}}={R^{n+1,m}}..................................Eq.(6)\)
Where
\(a_{w}^{{n+1,m}}= - \frac{{{k_w} \times \Delta t}}{{\rho \times \Delta X \times \delta {X_w}}},\) \(a_{e}^{{n+1,m}}= - \frac{{{k_e} \times \Delta t}}{{\rho \times \Delta X \times \delta {X_e}}},\)
\(a_{p}^{{n+1,m}}=\left( {C{{\left( T \right)}^{n+1,m}}+a_{w}^{{n+1,m}}+a_{e}^{{n+1,m}}} \right)\)
\({R^{n+1,m}}=\left( {h_{p}^{n} - h_{p}^{{n+1,m}}+C{{\left( T \right)}^{n+1,m}} \times T_{p}^{{n+1,m}}} \right)\)
The discretized Eq. (2–8) can be solved using an iterative correction scheme developed by Swaminathan and Voller [14], considering the PCM hysteresis in this case called: ICSH-TDMA algorithm. Noting that this method has been addressed by Al-Saadi and Zhai [11], and concluded that this scheme provides important advantages, is the flexibility of use by large steps, computational efficiency and stability in numerical predictions. Following the process of solving Eq. (6) is written as follows:
\({\left[ A \right]^{n+1,m}} \times {T^{n+1,m+1}}={R^{n+1,m}}...........................................Eq.(7)\)
In this numerical approach, after solving the equation for Tn+1, m+1 by using a direct solver based on the tri-diagonal matrix algorithm (TDMA), the enthalpy of the node is updated based on Eq. (4). At this moment, the enthalpy is known and therefore the temperature field is corrected to be consistent with the enthalpy temperature performance curve using the following relationship: [15]
\({T_p}=\left\{ \begin{gathered} \frac{{{h_p}}}{{{C_s}}}...........................................................,{h_p} \leqslant {C_s} \times \left( {{T_m} - \varepsilon } \right) \hfill \\ \frac{{{h_p}+\left[ {\frac{{{C_l} - {C_s}}}{2}+\frac{L}{{2 \times \varepsilon }}} \right]}}{{\left[ {\frac{{{C_l} - {C_s}}}{2}+\frac{L}{{2 \times \varepsilon }}} \right]}} \times \left( {{T_m} - \varepsilon } \right).............,{C_s} \times \left( {{T_m} - \varepsilon } \right)<{h_p}<{C_l} \times \left( {{T_m}+\varepsilon } \right)+L..................Eq.(8) \hfill \\ \frac{{{h_p} - \left( {{C_s} - {C_l}} \right) \times {T_m} - L}}{{{C_l}}}............................,{h_p} \geqslant {C_l} \times \left( {{T_m}+\varepsilon } \right)+L \hfill \\ \end{gathered} \right.\)
The iteration process continues in this prediction correction cycle until the convergence point is reached. It should be noted that the taking into account of the partial hysteresis of the phase change was adopted in the same way as that proposed by Bony and Citherlet [9]. the update of the parameters depends on the enthalpy in the algorithm, in the boxes framed in red as shown in Fig. 4, is done according to a new heat capacity graph as shown in Fig. 5. PCM phase change hysteresis was considered using a heat capacity curve in terms of two melting points, 22 ° C in the case of melting and 20 ° C in the case of solidification according to the proposed model by Feustel [16], Eq. (9).
\({C_p}\left( T \right)={C_{p,const}}+\frac{{{h_2} - {h_1}}}{2} \times \frac{{\frac{{2\beta }}{\tau }}}{{{{\cosh }^2}\left[ {\frac{{2\beta }}{\tau }\left( {T - {T_m}} \right)} \right]}}....................................................Eq.(9)\)
Where Cp = specific heat [kJ/kg K], T = Temperature [K], h = specific enthalpy [kJ/kg],
β = inclination [-], τ = width of the melting zone [K], Tm = melting temperature [K].
After solving the heat transfer equation and calculating the enthalpy, the algorithm checks the new enthalpy values, where if (hnew ≥ hold) the ICSH-TDMA program records the same values to do the next iteration, if (hnew<hold) the ICSH-TDMA algorithm corrects the enthalpy values according to the heat capacity curve (solidification curve). Before calculating the final temperature and storing it in the data.