**3. New Simulation Algorithm**

The hysteresis e ffect is usually shown in the form of two curves describing the course of the same phenomenon. The di fficulty associated with modeling the hysteresis is due to the need to introduce two di fferent functions describing the material enthalpy within the phase change range. The selection of the proper curve (Figure 1) must be based on the history of changes and the identification of the direction of the current changes. Such an algorithm would be appropriate if it could be assumed that the phase transitions occurred in the entire volume of the material.

**Figure 1.** Phase change hysteresis process represented as a function of temperature.

As outlined in the previous section, the real course of the phase changes is more complicated, and, in many cases, the process of melting/solidification is incomplete. A hysteresis algorithm should enable the direct transition from the melting curve to the freezing curve and vice versa. The concept of a 'state' was used by Zastawna [35] to create a new algorithm. Figure 2 shows the transition from the 'cooling curve state' to the 'heating curve state' defined as the occurrence of the third kind of state called the 'transition curve state'. The states are represented by the colors: blue—cooling curve; red—heating curve; purple-transition curve.

**Figure 2.** Intermediate transition curve between solidification (**a**) and melting (**b**) curves.

When the cooling phase of the completely melted material has begun (curve 0–1–2, Figure 2a) and before the phase change starts (point 2), the heat is released in a sensible way. After the phase change process has started, the energy is released in both latent and sensible ways (curve 2–3). At this point (point 3), the heating of the material begins again. Thus far, the enthalpy value follows the cooling curve state. The re-start of heating results in the termination of the cooling curve state and the start of the transition curve state (curve 3–4). For the transition curve, the enthalpy value is determined from either the heating curve or the cooling curve.

Rose et al. [36] suggested an instantaneous transition between the curves and a horizontal movement of the state until it reaches either the melting or solidifying curve. Bony et al. [37] assumed a sloped transition curve between the heating and cooling curves but a horizontal transition for the sub-cooling effect. Any change in the system temperature during the transition phase should obviously be associated with a change in enthalpy. The method of transition in the form of a horizontal line adopted in the article is a simplification.

While the hysteresis curves can be measured, the shape and slope of the transition curve are completely unknown, and experimental data are not available as they depend on many parameters, including the current entropy value. Curves describing the full transformation hysteresis are also a kind of idealization of this phenomenon and only describe it in an approximate way. Among others, heating and cooling curves are not usually parallel to each other, so any decision regarding the transition curve slope would always be approximate. Because the impact of simplification in the form of a horizontal transition line is small, and because of this, the modeling algorithm becomes significantly simpler, it was decided to use a horizontal transition curve in the model.

At point 4, the process of melting of the partially solidified material restarts. At this moment, the transition curve state ends, and the heating curve state will restart (curve 4–5 in Figure 2b). After completion of the phase change, the material again accumulates heat only in a sensible way (curve 5–1) while remaining in the heating curve state.

Equivalent specific heat of a PCM in case of the heating curve is calculated in the following way:

$$\mathcal{L}\_p(T) = \frac{H\_{\rm i,h}^j - H\_{\rm i,h}^{j-1}}{T\_i^j - T\_i^{j-1}} \tag{2}$$

where *Hj i*,*h* is the heating enthalpy of a discrete node '*i*' at a time-step '*j*' (and respectively '*j* − 1') [J/kg]. In the same way for the cooling curve:

$$\mathbf{c}\_p(T) = \frac{H\_{i,c}^j - H\_{i,c}^{j-1}}{T\_i^j - T\_i^{j-1}} \tag{3}$$

where *Hj i*,*<sup>c</sup>* is the cooling enthalpy of a discrete node '*i*' at a time-step '*j*' (and respectively '*j* − 1') [J/kg]. Forthetransitionfromcoolingtoheatingcurve3–4,Equationor,if*Hj* >*Hj* then:

$$\mathcal{L}\_P(T) = \frac{H\_{i,tr}^j - H\_{i,tr}^{j-1}}{T\_i^j - T\_i^{j-1}} \tag{4}$$

*i*,*tr* *i*,*h*

where *Hj i*,*tr* is the transition enthalpy, taken from transition curve and *Hj i*,*tr* = *Hj*−<sup>1</sup> *i*,*tr* If *Hj i*,*tr* ≤ *Hj i*,*h* then

$$c\_p(T) = \frac{H\_{\rm i,h}^j - H\_{\rm i,tr}^{j-1}}{T\_i^j - T\_i^{j-1}} \tag{5}$$

For the transition from heating to cooling curve 4–3, Equation or, if *Hj i*,*tr*< *Hj i*,*<sup>c</sup>*then:

$$x\_p(T) = \frac{H\_{i,tr}^j - H\_{i,tr}^{j-1}}{T\_i^j - T\_i^{j-1}} \tag{6}$$

and *Hj i*,*tr* = *Hj*−<sup>1</sup> *i*,*tr* If *Hj i*,*tr* ≥ *Hj i*,*<sup>c</sup>* then:

$$\mathbf{c}\_{\mathcal{P}}(T) = \frac{H\_{i,\mathcal{E}}^j - H\_{i,tr}^{j-1}}{T\_i^j - T\_i^{j-1}} \tag{7}$$

Appendix A shows the data that should be entered into the newly created PCM hysteresis algorithm. Each temperature corresponds to two enthalpy values separately for the heating and for the cooling curves.
