Reviewing Theoretical and Numerical Models for PCM-embedded Cementitious Composites

Abstract

Accumulating solar and/or environmental heat in walls of apartment buildings or houses is a way to level-out daily temperature differences and significantly cut back on energy demands. A possible way to achieve this goal is by developing advanced composites that consist of porous cementitious materials with embedded phase change materials (PCMs) that have the potential to accumulate or liberate heat energy during a chemical phase change from liquid to solid, or vice versa. This paper aims to report the current state of art on numerical and theoretical approaches available in the scientific literature for modelling the thermal behavior and heat accumulation/liberation of PCMs employed in cement-based composites. The work focuses on reviewing numerical tools for modelling phase change problems while emphasizing the so-called Stefan problem, or particularly, on the numerical techniques available for solving it. In this research field, it is the fixed grid method that is the most commonly and practically applied approach. After this, a discussion on the modelling procedures available for schematizing cementitious composites with embedded PCMs is reported.

1. Introduction

In the last century, the use of cementitious materials, like concretes and mortars, used as a base material for building and construction activities, has grown to vast quantities. The existing building stock, as well as the majority of newly established buildings and infrastructure, are nowadays largely made of reinforced and/or pre-stressed concrete. In the last decades, concrete innovations have been mainly focusing on improvements of the mechanical and durability properties, leading to systems with a high strength, high efficiency, and/or a high performance characteristics. Current endeavours are focussing on innovations that challenge the ability to reuse construction demolishing waste (CDW), saving energy demands while targeting a circular economy [1]. This trend has also been embraced by the concrete industry, which committed itself to reducing its carbon footprint dramatically and turning its “grey” image into a “greener” and more environmentally-friendly one [2]. Therefore, a significant effort still must be done in terms of enhancing the sector’s sustainability perspective by reducing its CO₂ footprint, widening its recyclability potential, and by cutting back on its energy demands [3].

Climate control and thermal comfort in modern buildings are the main issues that are of paramount importance for controlling and optimizing energy demands, total costs, and, as a spin-off, minimizing CO₂ emissions. However, implementing comfortability and, at the same time, sustainability in the new and existing building stock requires cutting edge innovations and efficient engineering systems, which are already the subject of current investigations [4]. Advanced technologies in the concrete industry make it possible to go even a step beyond these ambitions by turning cement-based materials into active thermal energy storing elements [5]. Thermal energy can be accumulated in a construction element using sensible, latent, and thermochemical-based materials. Among these, latent heat accumulating polymers, also known as so-called phase change materials (PCMs), are already a preferred solution in various applications. In this context, incorporating PCMs as a thermal energy storage material in construction elements has turned out to be a very efficient way to level out energy needs in modern buildings and/or other infrastructural applications. Due to their ability to absorb and liberate large quantities of latent heat at a constant temperature, PCMs have the potential to reduce the daily energy demand for heating and cooling, with a positive impact on global warming and climate change. Inclusion of PCMs in concrete leads, therefore, to a significant improvement of the material’s thermal performance. The efficiency to store thermal heat in a PCM-concrete makes it interesting to investigate the properties that allow for a wider use in construction and building applications [6,7].

Numerous different kinds of PCMs are being used in the literature. According to their chemical composition, they can be categorized as organic (O) compounds (paraffin and non-paraffin), inorganic (I) compounds (salt hydrate and metallic), and eutectic mixtures (O-O, I-I, I-O). PCMs can also be identified by their mode of phase transition: Liquid-gas, solid-gas, solid–liquid, and solid–solid. On the one hand, applications of both solid–gas and liquid–gas are limited in thermal storage capacity due to their large volume change, mainly associated with the gas state. On the other hand, the sensible heat capacity of solid-solid transitions is considerably lower than solid-liquid PCMs. For these reasons, in general, solid-liquid PCMs are preferred in thermal energy storage [8]. Other classifications of PCMs can be based on their thermal properties (e.g., melting temperature range, latent heat per unit volume, specific heat for sensible storages, thermal conductivity of their phases, etc.). For example, based on their temperature ranges, PCMs can be subdivided into low temperature PCMs (melting point below 10/15 °C) employed for cold applications, mid temperature PCMs (melting point between 15/80 °C) suitable for applications in building design, and finally, high temperature PCMs (melting point over 80 °C), normally used for industrial applications, like solar energy power plants [9].

A wide range of applications deal with the use of PCMs for enhancement of the thermal comfort of buildings and constructions. PCMs can be integrated in structural or non-structural elements for improvement of the thermal efficiency and sustainability in buildings in terms of cooling, heating, or both [10]. They can also be employed through passive or active systems for heating/cooling exchangers. In passive systems, the release of heating or cooling energy only depends on the environmental temperature fluctuations. These applications can easily be implemented in building frames, like walls, roofs, and/or floors, without any mechanical interference needed for changing or discharging heat. Contrarily, active systems do need a mechanical action to charge/discharge accumulated heat energy in PCMs [11].

Firstly, in Section 2, this state of the art document reports an overview of the available scientific literature on thermal energy storage and accumulations using PCMs. Section 3 reports a short discussion of the various experimental campaigns used for PCM-concretes, then, in Section 4, a discussion is addressed that focuses on a review of the various numerical tools used for modelling phase change problems, while emphasizing the so-called Stefan problem. Further details about the fixed grid method, adopted as the most common and practical approach for modelling the Stefan problem in PCMs, are described in Section 5. Then, Section 6 focuses on reviewing the available modelling procedures for cementitious composites enclosing PCMs based on different scale approaches (structural-, macro-, meso-, and microscale procedures). Finally, concluding remarks and an outline of future research steps to cover the observed research gap on numerical simulations of phase-change materials and processes in PCM-concretes will be reported in Section 7.

2. Review on Thermal Energy Storage Using PCM

Numerous articles on phase change materials have been reported in the literature in several fields of applications, with most of them of a scientific nature. Among several of them, the state of art of PCMs was explicitly addressed for emphasizing the use of thermal energy storage in various engineering and practical applications. These turned out to be very helpful for understanding the state of current research, and to determine the direction of future research steps for PCM developments.

This section reviews the available papers that use PCMs as a sensible latent storage system. Particularly,

Table 1. Latent heat storages with PCMs: General and buildings applications.

highlights those articles that explicitly address PCMs for thermal energy storages (TESs) usage. Until now, 23 works [9,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33] were found that deal with the use of PCMs as TES from a general perspective, such as industrial applications, buildings and construction elements, heating/cooling of water, solar energy plants, biomaterials and biomedical applications, electronic, automotive industry, food industry, etc. These articles focused on several aspects concerning the way TES was achieved by adopting PCMs, while accounting for different topics and/or focus areas, such as PCM types and classifications, sensible vs. latent TES in various amounts of PCMs, melting points, porous containers for PCMs, measurement techniques, experimental analysis, numerical tools and mathematical modelling, thermal conductivity enhancements, encapsulation methods, long term and cycle responses, cases studies, and applications.

Moreover, 37 articles [10,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68] dealt with PCMs for building applications. Among them, 23 referred to the available studies and research aimed at controlling the thermal comfort in buildings by means of PCMs [10,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55], seven focused on the employment of PCMs mainly in cooling applications [11,56,57,58,59,60,61], three were fully dedicated to the available theoretical models and numerical tools for modelling PCMs and PCMs in building applications [62,63,64], two were on microencapsulated (M)-PCM [65,66] in buildings, and, finally, the last two articles were on the thermal stability and thermal conductivity enhancement of PCMs [67,68]. The reviewed papers on PCMs in buildings, using latent thermal energy storage systems, covered several topics ranging from passive to active systems, heating and cooling applications for thermal comfort, experimental studies, small laboratory studies and full-size room tests, incorporation and impregnation techniques, PCM-concretes, sensible, latent and thermochemical-based PCMs for construction, enhancement solutions, micro and macro encapsulations, long term stability and cycle behavior, fire retardation in construction materials, the Stefan problem and numerical models, temperature-based heat capacity method, enthalpy formulations, fist law and second law based models, model collections, etc.

3. Experimental Investigations on PCM-Concretes

This section reports a short overview on the available experimental studies on PCM-concrete systems. The review is elaborated in a very broad way, thus, considering all applications with organic and inorganic PCMs. However, it may be worth mentioning that most PCM-concrete applications employ paraffin-based (organic) materials and only a few refer to other phase-change materials.

Numerous experimental-based research, available in the literature, addressed the hydro-thermo-chemo-mechanical properties of cementitious materials containing PCMs. Whenever PCMs are embedded in concrete used for passive heat storage in residential buildings, they are characterized by their melting temperature, which varies between 18 °С and 25 °С, which somehow corresponds to a standard temperature range for comfortable living [34]. PCMs employed in concrete can be of an organic (paraffin and non-paraffin) and/or inorganic (hydrated salts) nature [69]. The organic PCMs that are commonly employed for concrete are paraffin waxes, which are characterized by a very large latent thermal heat storage density, have little (or quasi absent) undercooling, low vapor pressure, are soundless, thermal and chemical stable, show a lack of phase separation, have a self-nucleating behavior, have adaptable phase change temperatures, are environmentally harmlessness, have no unpleasant odor, are non-toxic, and are economically attractive [70]. Despite this, one of the critical points that needs attention is the low thermal that conductivity paraffin waxes have, which reduces the pace at which heat can be stored and/or released during a melting and crystallization phase change operation [71]. Furthermore, another weak aspect is the density change that occurs during a heating/cooling cycle (i.e., which runs in both ways, i.e., through a melting/solidification as well as through a solid/liquid phase change), and which causes a change of its actual volume to a large extent [70].

Four typical methods for incorporating PCMs in concrete are distinguished: i.e., (i) immersion of porous concretes/mortars, (ii) direct mixing of micro-encapsulated (M-)PCMs, (iii) macro-encapsulation, and (iv) aggregate impregnation.

The first technique is based on the immersion of liquid PCM into the dried open pore space of a cementitious composite, where the effectiveness of it depends on the permeability/absorption capacity of the porous medium, the operating temperature, and the type of the employed PCM. These aspects may strongly affect the timespan needed to fully immerse the concrete. The method is thus suitable for concrete bricks, blocks, or wallboards, which are capable of absorbing large amounts of PCM; however, leakage problems, especially after various thermal cycling, may occur that may impact the environment and flammability of the final product, which implicitly represents the main drawbacks of this method. Experimental works investigating this technique were published by [72,73,74,75].

Direct mixing of microencapsulated PCMs in concrete and mortars are proposed in many articles. M-PCM are thin shell-based particles containing PCMs, and are proposed to be mixed in a concrete while avoiding leakage phenomena, and, at the same time, increasing the fire resistance of the microencapsulated PCM [76]. However, it may be worth mentioning that M-PCMs may lower the mechanical properties of a concrete, especially when particles are broken during the mixing and/or for high dosages [52]. M-PCMs in Portland cement concrete and geopolymer-concrete (GPC) were investigated in [77]. The results highlighted that on the one hand, the addition of microcapsules increased the porosity, resulting in a reduced compressive strength for both concrete and geopolymer systems. On the other hand, a lower thermal conductivity and a higher thermal energy storage was achieved due to the addition of M-PCM in both systems. Commercially available M-PCM (paraffin) “powders”, with a melting temperature of 23 °C, were used by [6] in concrete mixtures for analyzing its thermal performances. This work remarked that the thermal efficiency of PCMs in building applications was limited due to the very low thermal conductivity of the PCM-concretes. Mechanical, fracture, and thermal tests were conducted by [78] on concrete made with powder M-PCMs (phase change temperature of 23 ± 3 °C). This work highlighted the reduction in strength with increasing PCM dosages, however, this effect could be mitigated through the application of matrix strengthening, like, for example, silica fume additions. Contrarily, fracture properties, studied by various authors, were substantially less affected, especially for lower dosages of PCM. These results also highlighted the possibility to control thermal stresses in restrained concrete elements, whenever phase transitions occur during hydration at early-ages. The contribution in [79] deals with a new technique to measure the thermal conductivity and specific heat capacity of PCM-concrete bricks under phase change processes. Moreover, in [80], the thermal effect of micro-encapsulated PCMs mixed in concrete floors was investigated, while considering volumes that represented a scaled living space in dwellings. The behavior of self-compacting concrete containing micro-encapsulated PCMs was studied in [81]. This contribution also showed that with increasing amounts of PCM, a self-compacting PCM-concrete composite was obtained with a lower thermal conductivity and an increased heat storage capacity, both significantly contributing to the thermal performance and energy savings. As was also reported in other works, a significant loss in strength was observed whenever M-PCMs were added to a self-compacting concrete. A large number of M-PCM capsules were destroyed during the mixing process, as observed from a micro-structural analysis conducted by the authors [81]. The release of paraffin wax into the surrounding matrix was assumed to be responsible for the loss of strength. Real size concrete cubicles, for analyzing the improved thermal inertia as well as the inner thermal comfort, were investigated by [82] considering M-PCMs with a melting point of 26 °C. In recent years, additional mixing techniques have been proposed using expanded perlite particles [83]. These porous materials are capable of absorbing paraffin PCMs and can be used in direct mixing. It is remarkably that the initial principle of encapsulating waxes for novel thermal energy storage was developed for impregnating paraffin waxes into hydrophobically coated expanded perlite granules. In addition, binary mixtures consisting of capric, myristic, lauric, and palmitic acids were mixed as bio-based fatty acid PCMs in a solid form during concrete mixing by [84].

Macro-encapsulation deals with the integration of PCMs throughout “large” containers, such as tubes, spheres, panels, and/or other packages [31,41,52]. Among other works, macro encapsulated PCMs incorporated in concrete walls were discussed in [85]. In this regard, the effect of the position of macro-encapsulated PCM containers in concrete walls was evaluated, with special attention paid to the indoor temperature and humidity levels in the tested real size rooms. A thermal analysis of a concrete roof, having vertical cylindrical holes filled with PCMs, was realized and tested by [86]. Conventional and alveolar bricks with macro-encapsulated PCM were proposed by [87] as a Mediterranean construction under real conditions. The authors demonstrated a very good thermal behavior of these solutions in terms of energy savings and technical viability.

Aggregate impregnation is on PCMs embedded in porous aggregates. Thermal energy storage aggregate (TESA), e.g., natural, lightweights or recycled ones, is another possible technique to produce concrete and/or mortars with PCMs [88,89]. Thermo-physical properties and the mechanical behavior of cement-based mortars with PCMs incorporated in light weight aggregates (LWAs) were investigated by [90]. In this campaign, different dosages of paraffin waxes were included by immersion and by employing two alternative procedures to dry the LWA surfaces, i.e., compressed air or oven dried. The output of this work showed that high thermal inertia can be obtained whenever mortars are used that have LWAs filled with PCMs. In addition, it turned out that the thermal energy response could be enhanced when the PCM in LWA were filled up to a certain degree of its volume fraction. After that, the response could decrease. Fine lightweight aggregates with four different kinds of PCMs were used in PCM-mortars while aiming to reduce the freeze/thaw damage induced by various temperature cycles [91]. In this study, the PCMs were embedded in the aggregates by capillary suction during 24 h and at a constant temperature of 40 °C. The compressive strength, isothermal calorimetry, and thermal property were also analyzed. Experimental studies focusing on the benefits of LWAs used as internal reservoirs to store PCMs, for higher and lower temperature PCMs and for different freeze/thaw cycles, were performed by [92]. TESAs in concrete were also investigated by [89]. Image analysis and mercury porosimetry helped to understand the efficiency of these systems, as aggregates with a large pore connectivity and transport ability can absorb more PCMs. Particularly, it turned out that the inner geometry of the aggregates´ porous structure may have a significant effect on their ability to absorb PCMs.

4. Phase Change Materials: Stefan Problem

Phase change materials for building applications are characterized by solidification or melting when phase transitions take place, i.e., from liquid to solid or vice versa. In this particular situation, a heat-conduction problem according to the so-called Stefan problem can be applied.

The Stefan problem assumes the existence of two different domains, one representing the solid region (Ω_S) and another representing the liquid (Ω_L). The aforementioned domains are thus separated by a sharp surface front, namely a moving interface (Γ). Such a sharp surface front is characterized by a temperature field that equals the melting temperature, T_m (Figure 1).

The governing equations of the PCM heat conduction problem, under both solid–liquid (melting) and/or liquid–solid (solidification) transitions can be summarized as follows. Assumptions of the physical problem follow the classical approach of the Stefan problem, and considers it only as a conduction-based problem, a constant latent heat capacity (L), a fixed melting temperature for PCM (T_m), a zero-thickness interface of the sharp surface front (i.e., no mushy transition zone is considered), and the absence of nucleation and supercooling phenomena. According to these assumptions, the mathematical model can be formulated as follows.

Heat transfer problem in the solid (S) part:

$$\rho C_S\frac{\partial T_S}{\partial t}=\nabla .\left({\lambda }_S\nabla T_S\right)+{\dot{q}}_{v,S}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\forall x\in {\Omega }_S$$

(1)

Heat transfer problem in the liquid (L) part:

$$\rho C_L\frac{\partial T_L}{\partial t}=\nabla .\left({\lambda }_L\nabla T_L\right)+{\dot{q}}_{v,L}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\forall x\in {\Omega }_L$$

(2)

Heat balance at the sharp front interface (namely the Stefan condition):

$${\lambda }_S\nabla T_S\cdot n-{\lambda }_L\nabla T_L\cdot n=\rho L\frac{d\overline{\mathbf{x}}}{dt}\cdot n\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\forall \overline{\mathbf{x}}\in \Gamma$$

(3)

with $\rho$ as the density (considered constant when limited temperature variations occur, as is the case for PCM thermal problems, or, more generally, the density may depend on the position, X, and temperature field), ${\lambda }_S={\lambda }_S\left(T_S,x\right)$ and ${\lambda }_L={\lambda }_L\left(T_L,x\right)$ are the solid and liquid thermal conductivities, $C_S=C_S\left(T_S,x\right)$ and $C_L=C_L\left(T_L,x\right)$ are the heat capacities, $T_S$ and $T_L$ are the temperatures of the solid and liquid parts, while ${\dot{q}}_{v,S}$ and ${\dot{q}}_{v,L}$ are the source terms. Additionally, t represents the time variable, while ∇ and ∇ are the divergence and gradient tensorial operators, respectively. In Equation (3), $n=n(\overline{\mathbf{x}})$ is the unit normal vector on the moving interface, $\frac{d\overline{\mathbf{x}}}{dt}$ is the velocity of the interface, and L is the latent heat per unit mass.

Several authors have proposed analytical solutions for the above model, under the assumption of a simplified mono-dimensional situation for the heat conduction in the solid and liquid phases [93,94]. These proposals deal with simplified solutions for a one-dimensional model to determine the solid and liquid temperature field analytically as well as the solid–liquid interface during the solidification/melting processes (see a.o. the following works [95,96,97,98,99,100]).

More complex models can be handled by neglecting the aforementioned simplifications as were done in the original Stefan problem. In many solutions, the mushy zone (part of the materials’ interface in which the liquid and solid phases both may coexist) is not considered. Available proposals in the literature accounted for an explicit consideration of the mushy region. A computational fluid dynamics model was proposed by [101] for analysis of the phase change materials in building applications that consider three phase states, viz. mushy, solid, and liquid. The enthalpy approach outlined by [102,103] showed the ability to consider the possible occurrence of PCM mushy zones. These proposals were conceived in the general framework of convection/diffusion problems and showed compatibility such that they could be accounted for in fluid flow software. The same procedure, also known in the literature as the enthalpy-porosity approach, was followed recently by [104] for transient numerical simulations of melting/solidification processes of sodium nitrate (NaNO₃), used as a phase change material. In addition, the effect of internal air void formation in encapsulated PCMs were investigated by [105]. In this proposal, a combination of the volume to fluid method, following the enthalpy-porosity approach, were conducted to schematize all phases present in a capsule: i.e., solid, liquid, mushy part, and air.

Several authors pointed out the importance to also consider the conduction processes jointly with natural convection (mainly interesting for the motion of the liquid phases), which could significantly affect the phase change process. In this case, the governing equations for the convection-condition as a transient analysis of a melting process of PCMs follow the so-called Navier-Stokes equation, continuing equilibrium, and the energy equation [106,107]. However, it may be worth mentioning that if a thermal system contains a small fraction of PCMs (as is the case for cementitious composites with PCMs), only the conduction mode can be considered, which avoids the complexity of modeling the natural convection flow in the mushy and liquid phases [108,109,110].

The numerical tools needed for addressing the solution of the Stefan problem (and/or its extended formulations) can generally be subdivided into the following three different types:

• Fixed grid method: In this approach, a grid of spatial nodes used for discretizing the problem remains fixed during time, while the problem is studied with auxiliary constitutive formulations and state functions, which help to trace the phase change phenomena governed by solidification and/or melting. These models are mostly employed and several examples can be found in the scientific literature. Pioneer contributions in this fixed grid method can be found in [94,111,112,113].
• Deformed grid method: In this approach, the nodes forming the grid may move to follow the sharp fronts and interfaces occurring during melting and solidification. These models are thus able to solve the Stefan problem and to follow, by deforming the grid, i.e., the Stefan condition, as the solution evolves [114,115].
• Hybrid fixed/deformed method: These approaches couple solutions from discretized partial differential equations solved in a “ghost” fixed grid, with local front tracking schemes, which are based on deformed grid procedures [94,97,116].

The following section presents an overview of a fixed grid method, highlighting the governing equations and numerical solutions, and emphasizes the difference between them.

5. Fixed Grid Numerical Method

As already mentioned above, the fixed grid method is the simplest method to solve thermal problems involving phase changes. The basic equation of a heat conduction problem (the discussion is limited to conduction problems only, as they are applicable for cement-based composites in building physics problems) can be written as follows:

$$\frac{\partial H}{\partial t}=\nabla .\left(\lambda \nabla T\right)+{\dot{q}}_v\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\forall x\in \Omega$$

(4)

where H is the enthalpy of the system, λ = λ(T, x) is the thermal conductivity of the material, depending on the temperature, T, and position, X, and the heat capacities, while ${\dot{q}}_v$ is the source term.

The initial and boundary condition to solve Equation (4) are the following:

$$T\left(x,t=0\right)=T_0\left(x\right)\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\forall x\in \Omega$$

(5)

$$\begin{array}{ll}T\left(x,t\right)=T_D& \forall x\in {\Gamma }_T\\ \left(\lambda \nabla T\right)\cdot n=q& \forall x\in {\Gamma }_q\\ \left(\lambda \nabla T\right)\cdot n=h\left(T_{\infty }-T\right)& \forall x\in {\Gamma }_c\end{array}$$

(6)

${\Gamma }_T$, ${\Gamma }_q$ and ${\Gamma }_c$ are the boundaries of the domain $\Omega$, where the essential (Dirichlet) and natural boundary conditions are imposed; T_D is the specified temperature; q is the flux; h is the heat convection coefficient; and $T_{\infty }$ is the environmental temperature.

When the domain splits in two or three phases (i.e., liquid, solid, mushy), additional boundary conditions (Stefan conditions) need also to be considered. Particularly, in the case of an isothermal phase change, these additional conditions are represented by Equation (3) plus the following one:

$$T_S\left(x=\overline{\mathbf{x}}\right)=T_L\left(x=\overline{\mathbf{x}}\right)=T_m\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\forall \overline{\mathbf{x}}\in \Gamma$$

(7)

and analogously, in the case of a non-isothermal phase-change, which deals with the presence of a so-called mushy range, the Stefan conditions must be imposed, and can be extended to account for a third phase (mushy), where the heat balance is then imposed at two sharp front interfaces, i.e., a mushy-liquid surface and a mushy-solid one.

In the framework of the fixed grid method, several proposals can be found in the literature, always starting from Equation (4). From this, a distinction can be made mainly in three families: i.e., (i) the enthalpy-based method (EM), (ii) the apparent calorific capacity method (ACCM), and (iii) the heat source method (HSM). All these (sub-)methods are capable of determining a posteriori the moving front of melting or solidification. The following subsections describe a flavor of each of these approaches and highlight the available numerical techniques necessary for their solution.

5.1. Enthalpy-Based Method

The enthalpy-based method (EM), which solves a diffusion phase change problem in a general way, is based on a direct solution of the system of Equations (4)–(7).

In a general sense, the enthalpy of a system, H, in Equation (4) can defined as follows:

$$H\left(T\right)={\displaystyle \underset{T_{ref}}{\overset{T}{\int }}\rho C(T)dT}$$

(8)

where $\rho C(T)$ is the volumetric heat capacity and T_ref is a reference temperature.

In the case of a PCM, Equation (8) modifies into:

(i) PCM with an isothermal phase change:

$$H\left(T\right)=\{\begin{array}{cc}{\displaystyle \underset{T_{ref}}{\overset{T}{\int }}\rho C_S(T)dT}& \mathrm{when}T\le T_m(\mathrm{solid}\mathrm{phase})\\ {\displaystyle \underset{T_{ref}}{\overset{T_m}{\int }}\rho C_S(T)dT}+\rho L+{\displaystyle \underset{T_m}{\overset{T}{\int }}\rho C_L(T)dT}& \mathrm{when}T>T_m(\mathrm{liquid}\mathrm{phase})\end{array}$$

(9)

(ii) PCM with an non-isothermal phase change (presence of a mushy zone):

$$H\left(T\right)=\{\begin{array}{cc}{\displaystyle \underset{T_{ref}}{\overset{T}{\int }}\rho C_S(T)dT}& \mathrm{when}T\le T_{sol}(\mathrm{solid}\mathrm{phase})\\ {\displaystyle \underset{T_{ref}}{\overset{T_{sol}}{\int }}\rho C_S(T)dT}+\rho Lf(T)& \mathrm{when}{T}_{sol}<T\le T_{liq}(\mathrm{Mushy}\mathrm{phase})\\ {\displaystyle \underset{T_{ref}}{\overset{T_{sol}}{\int }}\rho C_S(T)dT}+\rho L+{\displaystyle \underset{T_m}{\overset{T}{\int }}\rho C_L(T)dT}& \mathrm{when}T>T_{liq}(\mathrm{liquid}\mathrm{phase})\end{array}$$

(10)

where T_sol and T_liq are the solid-mushy and mushy-liquid temperatures, respectively, while $\rho C_S(T)$ and $\rho C_L(T)$ are the volumetric heat capacity for the solid and liquid phases, respectively. The function, $f(T)$, in the mushy region is equal to 0 when T = T_sol and 1 when T = T_liq, allowing the definition of L (latent energy heat per unit mass), which represents the energy for a unit volume of solid to be completely transformed into a liquid.

In the case of constant heat capacities and linear expressions of the $f(T)$ function, Equations (9) and (10) reduces into:

(i) PCM with an isothermal phase change (Figure 2a):

$$H\left(T\right)=\{\begin{array}{cc}\rho C_{S}T& \mathrm{when}T\le T_m(\mathrm{solid}\mathrm{phase})\\ \rho \left(C_S{T}_m+L+C_L(T-T_m)\right)& \mathrm{when}T>T_m(\mathrm{liquid}\mathrm{phase})\end{array}$$

(11)

 (ii) PCM with a non-isothermal phase change (presence of a mushy zone, see Figure 2b):

$$H\left(T\right)=\{\begin{array}{cc}\rho C_{S}T& \mathrm{when}T\le T_{sol}(\mathrm{solid}\mathrm{phase})\\ \rho \left(C_S{T}_{sol}+L\frac{T-T_{sol}}{T_{liq}-T_{sol}}\right)& \mathrm{when}{T}_{sol}<T\le T_{liq}(\mathrm{Mushy}\mathrm{phase})\\ \rho \left(C_S{T}_{sol}+L+C_L(T-T_{liq})\right)& \mathrm{when}T>T_{liq}(\mathrm{liquid}\mathrm{phase})\end{array}$$

(12)

This approach is commonly employed for solving diffusion phase change problems, and was proposed for the first time by Eyres et al. [111]. An implicit enthalpy scheme for a one-dimensional phase change problem applied to the solution of a binary alloy solidification problem was discussed by Voller [117]. The author mainly proposed a method to solve a one dimensional (1D) problem with an isothermal phase change as outlined in Equations (9) and (11). After that, the model was extended to solve convection/diffusion phase change problems and was developed and described by Voller et al. [103,118]. A detailed discussion on the enthalpy method and its derivation of the apparent heat capacity as well as the source based method was given in [119]. In the framework of the enthalpy method, a finite element based formulation for solving coupled nonlinear heat conduction/diffusion processes, including phase changes, was developed by [120]. In addition, a fully implicit finite difference method was proposed by Costa et al. [121]. The method was validated for conduction plus convection under 1D conditions, while for a two-dimensional case, the same calculations were performed for a conduction case study only. Simple unidimensional computational conduction models for isothermal PCM, enclosed in single rectangular and cylindrical containers, were highlighted in [122]. Using the Galerkin finite element method, Bhattacharya et al. [123] solved phase change problems for pure and multicomponent phases, using the enthalpy method and paying special attention to the role of the mushy region. A Newton-type of iterative procedure was proposed by Nedjar [124] to solve stationary and transient problems in the field of both classical finite difference schemes and finite element procedures. Commercial paraffin waxes, stored in vertical and rectangular containers, were simulated with a 1D enthalpy-based model [125], employed for modelling a wall subjected to a constant heating-rate acting from one side, with an airflow stream on the other side. The enthalpy method was also used for solving the so-called Stefan problem through a Neumann-type boundary condition and a finite difference approximation by [126]. The difference between the enthalpy method and the effective heat capacity approach was outlined by Lamberg et al., who obtained several validations of the numerical tools using the commercial Femlab/Comsol program [127]. A slightly modified version of the enthalpy method was then proposed by Vyshak and Jilani [128]. In this proposal, the temperature was decoupled from the liquid fraction field, where this latter was updated, explicitly, through the currently known temperature field.

For building applications with PCM, the enthalpy method has been employed in a few examples. It can be mentioned that the study of Izquierdo-Barrientos et al. [129] deals explicitly with the transient heat transfer through different walls/layers. The authors developed a 1D transient conduction model for typical winter/summer temperature fluctuations. The commercially available computational fluid dynamics code FLUENT, which is based on the enthalpy method, was employed by [130,131] to conduct numerical studies for the thermal performance of building walls and roofs with PCMs for residential applications. Finally, recently, an interesting numerical study based on the enthalpy method was employed in the medicine field by Kumar et al. [132]. Particularly, this work reported a numerical investigation of a phase change heat transfer process in a lung undergoing cryosurgery.

5.2. Apparent Calorific Capacity Method (ACCM)

The apparent calorific capacity method (ACCM) considers the effect of enthalpy and its evolution in time by considering an apparent (or sometime called effective) heat capacity during thermal phase changes.

The approach is based on the following relationship:

$$\frac{\partial H}{\partial t}=\frac{\partial H}{\partial T}\frac{\partial T}{\partial t}$$

(13)

where $\frac{\partial H}{\partial T}=\rho C_{eff}\left(T\right)$ represents the temperature-dependent apparent (effective) heat capacity [133].

According to this, Equation (4) modifies into the classical heat equation as follows:

$$\rho C_{eff}\left(T\right)\frac{dT}{dt}=\nabla .\left(\lambda \nabla T\right)+{\dot{q}}_v\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\forall x\in \Omega$$

(14)

The method has been well used in the literature since, in this approach, the temperature is the only variable that needs to be solved, discretized, and simulated. Actually, Equation (14) shows a problem that mainly deals with a classical heat problem. The key novelty of this method is to achieve a proper representation of the apparent capacity, C_eff.

Table 2. ACCM-based proposals in the framework of fixed grid methods.

Authors	Apparent Heat Capacity—Ceff
Lamberg et al. (2004) [127] Bonacina et al. (1973) [134] Idelsohn et al. (1994) [135] Voller (1997) [136] Pasupathy et al. (2008) [137,138] Zang et al. (2008) [139] Fang and Medina (2009) [140] Thiele et al. (2015) [141] Šavija and Schlangen (2016) [142]	$C_{eff}=\{\begin{array}{l}{C}_S\\ \frac{L}{2\Delta T}+\frac{C_S+C_L}{2}\\ C_L\end{array}\begin{array}{c}T\le T_m-\Delta T\\ T_m-\Delta T<T\le T_m+\Delta T\\ T>T_m+\Delta T\end{array}$
Hu and Argyropoulos (1996) [143]	$C_{eff}=\{\begin{array}{l}{C}_S\\ C_{in}\\ C_L\end{array}\begin{array}{c}T\le T_{sol}\\ T_{sol}<T\le T_{liq}\\ T>T_{liq}\end{array}$   $C_{in}=\frac{{\displaystyle \underset{T_{sol}}{\overset{T_{liq}}{\int }}C(T)dT+C_s{T}_{sol}}}{\left(T_{liq}-T_{sol}\right)}$
Samarskii and Vabishchevivh (1995) [144] Heim and Clarke (2004) [145] Heim (2010) [146]	$\frac{\partial H\left(T\right)}{\partial t}=\rho C_{eff}(T)\frac{\partial T}{\partial \upsilon }\frac{\partial \upsilon }{\partial t}$   $\upsilon ={\displaystyle \underset{C_{sol}}{\overset{C_{liq}}{\int }}{C}_{eff}(T)dT}$
Zukowski (2007) [147] Mankel et al. (2019) [148]	$\begin{array}{l}{C}_{eff}=C_{p,\mathrm{P}}=f\left(T_{\mathrm{P}}\right)\\ \mathrm{tabular}\mathrm{data}\mathrm{and}\mathrm{interpolated}\mathrm{with}\mathrm{cubic}\mathrm{spline}\mathrm{function}\mathrm{method}\end{array}$ tabular data from DSC tests
Evola et al. (2013) [149]	$\begin{array}{l}{C}_{eff}=\{\begin{array}{l}1200+18800\cdot e^{-\left(\frac{Tp-T}{1.5}\right)}\\ 1300+18700\cdot e^{-4{\left(Tp-T\right)}^2}\end{array}\begin{array}{c}T\le T_p\\ T>T_p\end{array}\\ \mathrm{equivalent}\mathrm{heat}\mathrm{capacity}\mathrm{of}\mathrm{PCM}\mathrm{wallboards}\end{array}$
Thiele et al. (2015) [141,150] (for composite applications: i.e., cementitious ones)	${\left(\rho C\right)}_{eff}\left(T\right)=\{\begin{array}{l}{\left(\rho C\right)}_{eff,S}\\ {\left(\rho C\right)}_{eff,S}+{\varphi }_c\frac{L}{2\Delta T}\\ {\left(\rho C\right)}_{eff,L}\end{array}\begin{array}{c}T\le T_m-\Delta T\\ T_m-\Delta T<T\le T_m+\Delta T\\ T>T_m+\Delta T\end{array}$   $\begin{array}{l}\underset{composite}{{\left(\rho C\right)}_{eff,S}}=\underset{coreoftheMPCM}{{\varphi }_c{\left(\rho C\right)}_{c,S}}+\underset{shelloftheMPCM}{{\varphi }_s{\left(\rho C\right)}_s}+\underset{concretematrix}{\left(1-{\varphi }_c-{\varphi }_s\right){\left(\rho C\right)}_m}\\ \begin{array}{cc}\begin{array}{l}{\varphi }_c\mathrm{volume}\mathrm{fraction}\mathrm{of}\mathrm{MPCM}\mathrm{core}\\ {\varphi }_s\mathrm{volume}\mathrm{fraction}\mathrm{of}\mathrm{MPCM}\mathrm{shell}\\ {\varphi }_m\mathrm{volume}\mathrm{fraction}\mathrm{of}\mathrm{concrete}\mathrm{matrix}\end{array}& {\varphi }_c+{\varphi }_s+{\varphi }_m=1\end{array}\\ \underset{}{{\left(\rho C\right)}_{eff,L}}=\underset{}{{\left(\rho C\right)}_{eff,S}}\end{array}$
Tittelein et al. (2015) [151] Kheradmand et al. (2016) [152] (for composite applications: i.e., cementitious ones)	$\begin{array}{l}{C}_{eff}=X_{mort}\cdot C_{mort}+X_{PCM}\cdot C_{PCM}\\ \mathrm{being}{X}_{mort}\mathrm{and}{X}_{PCM}\mathrm{mass}\mathrm{fractions}\end{array}$   $\begin{array}{l}{C}_{mort}=\mathrm{constant}\\ C_{PCM}=f\left(T\right)\end{array}$
Yao et al. (2018) [153]	$C_{eff}=\{\begin{array}{l}\frac{C_S+C_L}{2}\\ \frac{L}{2\Delta T_1}+\frac{C_S+C_L}{2}\\ \frac{C_S+C_L}{2}\end{array}\begin{array}{c}T\le T_m-\Delta T_1\\ T_m-\Delta T_1<T\le T_m+\Delta T_1\\ T>T_m+\Delta T_1\end{array}$ melting $C_{eff}=\{\begin{array}{l}\frac{C_S+C_L}{2}\\ -\frac{L_2}{2\Delta T_2}+\frac{C_S+C_L}{2}\\ \frac{C_S+C_L}{2}\end{array}\begin{array}{c}T\le T_s-\Delta T_2\\ T_s-\Delta T_2<T\le T_s+\Delta T_2\\ T>T_s+\Delta T_2\end{array}$ solidification
Young et al. (2018) [154]	$\begin{array}{l}\frac{{\left(\rho C\right)}_{eff}\left(T\right)}{{\left(\rho C\right)}_{eff,s}}=\{\begin{array}{l}1\\ 1+{\varphi }_{c+s}\frac{L^{*}}{\Delta T}\\ 1\end{array}\begin{array}{c}T\le T_m-\frac{\Delta T}{2}\\ T_m-\frac{\Delta T}{2}<T\le T_m+\frac{\Delta T}{2}\\ T>T_m+\frac{\Delta T}{2}\end{array}\\ L^{*}=\mathrm{dimensionless}\mathrm{latent}\mathrm{heat}\mathrm{of}\mathrm{fusion}\end{array}$
Singh and Bhat (2018) [155]	$C_{eff}=\{\begin{array}{l}{C}_S+C_{p,\max }\cdot e^{\left(\frac{-{\left(T_m-T\right)}^2}{2W^2}\right)}\\ C_L+C_{p,\max }\cdot e^{\left(\frac{-{\left(T_m-T\right)}^2}{2W^2}\right)}\end{array}\begin{array}{c}T\le T_m\\ \\ T>T_m\end{array}$
Han et al. (2018) [156]	$\begin{array}{l}\underset{composite}{{\left(\rho C\right)}_{eff}}=\underset{phase1}{{\varphi }_{ph1}(T){\rho }_{ph1}{C}_{ph1}}+\underset{phase2}{\left(1-{\varphi }_{ph2}(T)\right){\rho }_{ph2}{C}_{ph2}}+\underset{latentpart}{\rho L\frac{\partial {\alpha }_m}{\partial T}}\\ {\alpha }_m=\mathrm{smooth}\mathrm{transition}\mathrm{function}\end{array}$

outlines various proposals for describing the apparent calorific capacity for a phase change process as available in the scientific literature.

5.3. Heat Source Method (HSM)

The heat source method (HSM) separates the enthalpy into two parts, viz. a sensible and a latent one. The latent term is thus used as a source term in the classical heat equation as follows:

$$\frac{\partial H}{\partial t}=\frac{\partial H}{\partial T}\frac{\partial T}{\partial t}=\left[\rho C(T)+\rho L\frac{\partial f(T)}{\partial T}\right]\frac{\partial T}{\partial t}$$

(15)

Combining Equations (4) and (15), the following heat equation can be obtained:

$$\rho C\left(T\right)\frac{dT}{dt}=\nabla .\left(\lambda \nabla T\right)+{\dot{q}}_v-\rho L\frac{df(T)}{dt}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\forall x\in \Omega$$

(16)

where $-\rho L\frac{df(T)}{dt}=S$ represents the so-called source term due to the latent part of the enthalpy.

This method was originally proposed by Eyres et al. [111] based on the idea that heat enthalpy, H, is function of the temperature and volume fraction of the material that is affected by the phase change process. According to this hypothesis, the author reached an expression similar to Equations (15) and (16), envisioning the pioneering contribution of this proposal. An alternative proposal was developed by Voller (1995) [157], where, based on an implicit finite-difference discretization, the enthalpy formulation was separated into a sensible and latent heat term. The resulting set of non-linear equations was based on isolation of the source latent term while the method offered the advantage that only the temperature field needs to be solved in the resulting iterative heat equation scheme. Following this proposal, a new coupled HSM-based scheme was developed by Swaminathan and Voller (1997) [158] aimed at simulating the solidification behavior of multicomponent alloys for a wide range of local scale behaviors. Latent heat storage systems were also designed by Costa et al. (1998) [121], using an enthalpy formation solved through a fully implicit finite difference method and the HSM. The proposal of these authors deals with a non-linear system of equations with the only isolated variable being the enthalpy (so, not the temperature as usually adopted in these thermal-based problems). The solution of a temperature-based transient heat conduction problem involving phase changes was solved by Fachinotti et al. (1999) [159] in the framework of a finite element model and using the HSM approach. Brusche et al. (2006) [160] used the HSM for modelling phase-change layers in optical rewritable recordable media, while prior to that, a mathematical model was developed by Rostamizadeh et al. (2012) [161], showing an enthalpy formulation with the HSM technique and the effect of PCM thickness in energy storage systems. All the discussed proposals are summarized in

Table 3. HSC proposals in the framework of fixed grid methods.

Authors	Heat Capacity—C	Source Term—S
Voller (1995) [157]	$\rho C$	$-\rho \frac{\partial \Delta H}{dt}\Delta H=\{\begin{array}{l}0\\ L/2\epsilon \left(T+\epsilon \right)\\ L\end{array}\begin{array}{c}-\epsilon >T\\ -\epsilon \le T\le \epsilon \\ T>\epsilon \end{array}$
Swaminathan and Voller (1997) [158]	$\begin{array}{l}\rho C=(1-g){\rho }_S{C}_S+g{\rho }_L{C}_L\\ g=\mathrm{liquid}\mathrm{volumetric}\mathrm{fraction}\end{array}$	$\begin{array}{l}-{\rho }_S\Delta H\frac{\partial g}{dt}\mathrm{with}g=g\left(T\right)\\ \Delta H=\mathrm{constant}\end{array}$
Costa et al. (1998) [121] Athukorallage et al. (2018) [162]	${\rho }_S{C}_S\mathrm{and}{\rho }_L{C}_L$	$\begin{array}{l}-{\rho }_{L}L\frac{\partial f}{dt}f=\{\begin{array}{l}0\\ 1\end{array}\begin{array}{c}T\le T_m\\ T>T_m\end{array}\\ 0\le f\le 1\mathrm{in}\mathrm{case}\mathrm{of}\mathrm{Mushy}\mathrm{region}(T_{sol},T_{liq})\end{array}$
Fachinotti et al. (1999) [159]	$\rho C$	$-\rho L\frac{\partial f}{dt}f=\{\begin{array}{l}0\\ 0<f^{*}\left(T\right)\le 1\\ 1\end{array}\begin{array}{c}T\le T_{sol}\\ T_{sol}<T\le T_{liq}^{}\\ T>T_{liq}\end{array}$
Brusche et al. (2006) [160] Rostamizadeh et al. (2012) [161]	$\rho C$	$-\rho L\frac{\partial f}{dt}f=\{\begin{array}{l}0\\ 1\end{array}\begin{array}{c}T\le T_m\\ T>T_m\end{array}$

.

5.4. Concluding Remarks on Fixed Grid Method

The above described methods are employed to simulate phase change processes governed by solidification and melting phenomena. As outlined before, in the enthalpy method (EM), the enthalpy variable is the decisive parameter that accounts for separating the sensible and latent heat parts for both isothermal phase changes and for those cases where mushy regions are present. Contrarily, the apparent calorific capacity method (ACCM) is a simpler method where the sensible and latent heat energy parts are lumped. Finally, the heat source method (HSM) treats the latent heat as a source term to obtain a temperature-based differential equation to be solved.

The EM and HSM are the only choices to be followed in the case of isothermal phase changes and/or for those cases where a phase change occurs in a small temperature range. On the other hand, ACCM is very easy to program being quite similar to the classical heat (conduction) problem with the only difference that the heat capacity is now temperature-dependent. ACCM is suitable for gradual phase changes, while it is not useful for isothermal phase changes (a jump in the H-T curve cannot be well captured by an appropriate description of C_eff).

In the framework of a finite element procedure, the preferred methods to be used are the ACCM and the HSM, due to the fact that the only dependent variable is temperature, which characterizes the problem. Even though, since the HSM can be employed for gradual, sharp, and isothermal phase processes, this method is the best choice for FE simulations and PCM analysis.

6. Cementitious Porous Materials with PCM: Numerical Models and Scales

Nowadays, many theoretical and numerical models are available for simulating heat storage phenomena in porous composites, like concrete, mortar, gypsum, etc. This section reports the available models developed in the last two decades and categorizes them by means of the different scales of observation.

Concrete and other cementitious materials are multiphase (composite) materials. For this reason, it can be considered and modelled as a homogeneous continuum at the macroscale and even structural scale, while at lower levels (meso-, micro-, and nano-scales), it can be considered as an inhomogeneous multiphase material.

6.1. Structural-Scale Models

These models allow capturing of the essence of heat storage phenomena at the structural (building physics) scale level. The major objective of these formulations is to build up simplified methods, sometime practice-oriented, for predicting the temperature fluctuations and thermal comfort in certain construction.

Numerical tools for analyzing multi-zone solar buildings, having PCM-impregnated gypsum plasterboard, were proposed by Heim et al. [145,146,163]. These authors used the commercial ESP-r program (numerical code for building physics performance simulation) to incorporate the effects of phase change materials. Further case studies concerning whole building thermal performances dealing with wallboards containing PCMs and micro-encapsulated PCMs in thin wallboards were investigated by Evola et al. in [149] and [164], respectively.

Moreover, building physics simulations using available commercial codes for PCMs in building applications are those related to the use of TRNSYS (Transient System Simulation Tool) [165,166,167], the International Building Physics Toolbox (libraries and models for Matlab/Simulink) [168,169,170,171], and other commercial programs, such as EnergyPlus, RadCool (a Web-enabled Simulation Tool for Radiative Cooling), CoDyBa, BSim, PCM Express and Wufi [172].

6.2. Macro-Scale Models

In the field of computational mechanics, macro-scale models for analyzing heat transfer phenomena in cement-based materials with PCMs are most popular and most frequently adopted. These models are based on the assumption that the schematized material acts as a continuum and homogenous medium. The most classical examples in this field are those related to the use of the aforementioned EM, ACCM, and HSM techniques [173] for macroscopic analysis of cementitious composites. Examples are from Pasupathy and Velraj [137,138], who analyzed the thermal performance of building roofs with PCMs. The authors employed the ACCM in a macroscale model where the concrete slab, the roof top, and PCMs were considered as homogeneous and isotropic elements. To study the influence of a multi-layer wall, including PCMs, a macroscale transient heat transfer model was developed by [129]. In this study, the authors focused on the effects of different boundary conditions, wall configurations and orientations, ambient conditions, and phase transition temperatures of the PCM. Then, a macroscopic enthalpy-based model was analyzed by [174] with the aim to evaluate the related errors occurring when misinterpreted enthalpy functions are employed in thermal simulations. Furthermore, a four-layer roof containing PCM was approached macroscopically by [131] using EM in a Fluent software. In this respect, the recent works based on the fixed grid methods can be cited, i.e., Kong et el. [130], Tittelein et al. [151], and Kheradmand et al. [152].

A one-dimensional analytical model was proposed by Mirzaei and Haghighat [101], where the concept of resistances and capacitances with computational fluid dynamics ingredients were combined for analyzing multi-layer cementitious walls with PCMs. Similarly, a simplified dynamic model for concrete brick-walls with PCMs, based on the aforementioned resistances/capacitances, and combined with genetic algorithm-based procedures was developed by Zhu et al. [175,176]. In the field of finite element approaches, heat transfer processes of gypsum boards containing fatty acid PCMs was studied by [177] using the Ansys FEM code. Then, an FEM-based approach considering the ACCM was followed by Thiele et al. [150], who proposed a technique for estimating the specific heat capacity for a PCM-composite.

Several homogenization techniques were then proposed for estimating the macroscopic heat properties of cement-based composites. An interesting proposal was given by [178] who considered the PCM-composite to be homogeneous and isotropic, and outlined a homogenization technique that explicitly considered the superficial capsule area per unit mortar volume (expressed as a function of the volume fraction of PCM to the surrounding mortar and to the capsule diameter), employed for evaluating the composite thermal properties of the material. PCM-mortar composites were schematized as a homogeneous medium by [179], hence considering the average and isotropic properties. These authors proposed a mixed rule for accounting the mass enthalpy of a composite material (i.e., cement mortar plus the PCM). Furthermore, a macroscale model for simulating concrete slabs cast with microencapsulated PCMs was proposed by [180], where the thermo-physical properties of the PCM-concretes were estimated through an application of mixture theory, while an analytical expression was introduced for the effective apparent calorific capacity.

Some recent works also accounted for the estimation of PCM effects on the hydration process, and in particular on the mitigation of the reaction heat resulting as a partial reaction of the cement grains with water. A theoretical model of a heat-flow and hydration model for concrete hardening, including paraffin, was presented by Pei et al. [181]. Furthermore, a proposed model for early-age PCM-concretes by Arora et al. [182] demonstrated significant reductions in the peak hydration temperature when PCMs were considered as the partial replacement of the cement paste or fine aggregates. The authors also proposed a formulation that represents the hydration-dependency of the PCM thermal conductivity and heat capacity. A study to examine the benefits related to the addition of microencapsulated PCMs in concrete pavement sections to mitigate the temperature rise during early-age cement hydration was recently proposed by Young et al. [183].

Equivalent thermal properties that can be used in the aforementioned macro-scale models and also the structural ones have been subject of research as well, with some works available in the literature. For example, a method based on a sequential quadric programming procedure to determine the specific calorific capacity and thermal conductivity of PCM-concrete bricks was proposed by [79]. A general discussion with many examples was done by Felske [184], addressing an evaluation of the effective thermal conductivity when composite spheres were merged in continuous media. The proposal by Felske [184] has been recently extended and verified by Thiele et al. [185] and Ricklefs et al. [186] for predicting the effective thermal conductivity of composite materials with PCM capsules. The Mori-Tanaka principles were then used by [187] for obtaining the macroscopic thermal conductivity of composites considering the effect of the interfacial particle resistance and their size distribution. The same approach was also followed by Stránsky et al. [188] to analyze the effective thermal conductivity of matrix-inclusion composites, and where some aspects were investigated, such as the orientation of particles and the perfect vs. imperfect interface. Furthermore, multi-step Mori-Tanaka homogenization techniques were employed by Aguayo et al. [189] to predict the thermal conductivity of the PCM-composites. Finally, a new hybrid technique was developed by Mandilara et al. [190] for estimating the effective heat capacity to be used in numerical simulations that account for PCM-concrete composites.

6.3. Meso-Scale Models

A better understanding of PCM mechanics subjected to external thermal fluctuations can be obtained by considering the mesoscale nature of the PCM composite materials. The mesoscale of concrete can be idealized by considering different phases, which together constitute such a composite. Thereby, the interaction among the different phases of the composite (i.e., matrix, coarse aggregates, PCMs, and interfaces between them) is explicitly considered in these approaches.

Only a few approaches are currently available in the literature for modeling meso-scale systems accounting for phase changes in concrete. In this field, a very simplified approach was proposed by Zhang et al. [139], where a 2D discretized model of a composite PCM-concrete was reported. The authors considered the PCM particles as circles with a radius of 1 mm and were uniformly distributed over the 2D geometry. No interface and thermal resistance were adopted. Furthermore, a heterogeneous meso-scale proposal based on a unit cell containing core-shell particles arranged in a periodic face-centered cubic packing was proposed by Thiele et al. [141]. This research was employed for validating a homogenization technique for continuous-based simulations compared with the same results obtained from a heterogeneous mesoscopic representation. Finite element simulations carried out on cementitious microstructures were proposed by Agauyo et al. [191] to investigate the influence of interface properties and PCM inclusions. The model explicitly considered the meso-scale with PCM cores and PCM shells of the micro-encapsulated PCM particles, the PCM-mortar interfaces, and the surrounding mortar matrix. A multi-scale model based on the analysis of a meso-scale RVE was also proposed by Mohaine et al. [192]. Particularly, the authors generated mesoscopic scale specimens based on a digital concrete model for generating random structures containing heterogeneities, like coarse aggregates and PCMs. The homogenization multiscale procedure, based on the finite element method, enabled the authors to obtain the thermo-mechanical parameters of cement paste and the concrete with PCMs. Furthermore, a two scale numerical model was proposed by Šavija and Schlangen [142], where at the meso-scale, the PCM microcapsules were explicitly simulated as embedded inclusions in the concrete matrix based on a discrete (lattice) model, while at the macro-scale, a commercial FE package was used to calculate the stress distributions from the PCM-affected temperatures in a wall-on-slab system. A meso-scale approach for determining the thermal deformation properties of microencapsulated PCMs embedded in cementitious composites was proposed by Young et al. [193], where various cell configurations with many different particle size distributions and packing arrangements were considered. Finally, the authors proposed a multiscale and thermodynamically consistency framework for solving thermo-mechanical problems with phase change phenomena [194].

7. Concluding Remarks

Accumulating solar and/or environmental heat in innovative materials for construction and building applications represents a way to level-out daily temperature differences and to significantly cut back on energy demands. A way to achieve this goal is to develop advanced cementitious composites that are characterized by a porous microstructure, which have the ability to accommodate phase change materials (PCMs) with the potential to store/relief energy during their phase changes from solid to liquid or vice versa.

From the current literature review, the following concluding remarks can be drawn:

• A large number of studies on various PCM-related topics are available in the literature. Some of them referring to specific topics from a material point of view (i.e., micro- and macro-encapsulation of PCM, thermal conductivity enhancements, thermal stability), while others refer to the possible field of application of PCMs, like thermal energy storages, high temperature applications, solar water heaters, cold applications, and building accumulation solutions;
• numerical solutions for analyzing the so-called Stefan problem in phase change materials were reviewed and discussed. The literature on different solution approaches was presented discussed;
• fixed grid method, representing the simplest approach for solving thermal problems involving PCMs, were reviewed and the basic equations behind this approach were outlined and discussed;
• in the framework of the fixed grid method, three different approaches were recognized, namely the “enthalpy-based method (EM)”, the “apparent calorific capacity method (ACCM)”, and the “heat source method (HSM)”. All these methods have in common that they are capable of determining a posteriori the movement of the melting/solidification front according to the Stefan problem. The different methods and corresponding models associated with the EM, ACCM, and HSM approach were briefly described and discussed; and
• the available models for PCM-concrete composites, focusing on the different length scales, i.e., micro, meso, macro, and possible multiscale approaches, are being reviewed and discussed. Many articles on PCM-concretes followed a macroscopic and building scale approach, while only few contributions are available in the literature that deal with the meso- to micro-scale level.

As a final comment, it may be worth mentioning that although a significant research effort has already been contributed to the field of numerical modelling for heat transfer processes with PCM accumulations, further efforts in the field of cementitious composites embedding PCMs are certainly embraced. Unambiguous knowledge on how porous microstructures, like cementitious composites or aggregates, should be built up and/or what are the critical demands that allow empty pores to serve as “closed” encapsulation cavities is still lacking. A holistic approach, combining a micro- to mesoscale poro-analysis model with a multiscale/multiphysics approach, along with a microstructural response, moisture diffusivity, phase change, and thermal analysis, has not been published in the literature so far.

Recent developments on numerical modelling of composite systems, including micro- and mesoscale observations, for considering local effects, like, for example, pores, air bubbles, PCMs, aggregates, interfaces, etc., are currently ongoing at the WiB Institute of TU-Darmstadt. Advantages of these lower-scale simulations is that they allow optimization of the design of efficient PCM-cement systems, thus avoiding the performance of time-consuming experimental tests. On the other hand, multiscale modelling applications are also under development for simulating at several scales the energy storage capacity of PCM-concretes. There, a multiscale framework, consisting of coupled mechanisms with thermal/mechanical and phase transformation phenomena, was considered.