Thermomechanical Performance Analysis of Novel Cement-Based Building Envelopes with Enhanced Passive Insulation Properties

The design of new insulating envelopes is a direct route towards energy efficient buildings. The combinations of novel materials, such as phase-change (PCM), and advanced manufacturing techniques, such as additive manufacturing, may harness important changes in the designing of building envelopes. In this work we propose a novel methodology for the design of cement-based building envelopes. Namely, we combined the use of a multiscale, multiphysical simulation framework with advanced synthesis techniques, such as the use of phase-change materials and additive manufacturing for the design of concrete envelopes with enhanced insulation properties. At the material scale, microencapsulated PCMs are added to a cementitious matrix to increase heat storage. Next, at the component level, we create novel designs for the blocks, here defined as HEXCEM, by means of additive manufacturing. The material and component design process is strongly supported on heat transfer simulations with the use of the finite element method. Effective thermal properties of the mixes can be obtained and subsequently used in macroscale simulations to account for the effect of the volume fraction of PCMs. From the experimental and numerical tests, we report an increase in the the thermal inertia, which results in thermal comfort indoors.


Introduction
The building sector is responsible for 40% of the primary energy requirements in the EU and 36% of greenhouse gas emissions in Europe [1,2]. Therefore, there is a global need for improving energy management in buildings, which can be approached using different strategies [3] dealing with generation, energy storage, materials, or control systems that not only are more sustainable and environmentally-friendly, but also economical. One solution for providing thermal comfort while reducing the energy demand is the use of new insulating building envelopes [4,5].
Active/passive thermal energy management systems [6,7] based on renewable sources (e.g., solar radiation) are gaining much attention within the nearly zero energy buildings (NZEBs) paradigm [8]. While active systems are generally more difficult to implement in a building, passive systems can be easily integrated in building envelopes [9], with neither aesthetic nor functional drawbacks [7].
Concrete is one the most widely used materials worldwide [10,11]. It is mainly used in buildings for structural elements (i.e., columns and beams), foundations, slabs, and envelopes. Although concrete has caught attention for its mechanical properties, it can also be designed with targeted thermal insulation properties. For instance, higher internal devices, such as thermal batteries [46], geothermal dissipation devices [47], freeze-thaw damage-resistant devices [48], or early-age cracking mitigation devices [16].

Materials and Methods
In this section we present the multiscale methodology followed for the design and validation of the building envelopes. In the first place, we designed a cementitious matrix with microencapsulated PCM at different volume fractions. These mixes were then characterized numerically using the FEM. The resulting material properties were subsequently upscaled. Then, we carried out the component design by considering different cell sizes. We used the FEM at the component scale for two purposes: (i) to validate the designs and (ii) to develop a multiscale simulation framework for new building envelopes.

Material Design
The material designed for the building envelopes consists of a CEM I 42.5 cementitious matrix (ENCI, The Netherlands) and microencapsulated PCM (Encapsys LLC, United States) with different volume fractions: V f = 0%, 10%, and 30%. The mix had a water-to-cement ratio of w/c = 0.45 and PCM capsules reach a random distribution within the cement paste, as shown in Figure 1. The PCMs used are organic paraffin encapsulated by a melamine-formaldehyde shell. The same PCMs have been used by the authors in reference [49]. Further details on the PCM characterization can be found in references [49,50]. Table 1 shows the composition of the three mixes.  [49]). PCMs from the same batch were used in the current work.
Once ready, the mixes were poured into the molds of the bricks and hardened for 2 days. Then, these were demolded and placed in a cure room until 21 or 28 days. According to [18], the thermal conductivity of cement paste decreases during the first 7 days of curing, but remains practically constant afterwards. The PCMs used are organic paraffin encapsulated by a melamine-formaldehyde shell, as in reference [16]. The heat of fusion, h f , is 143.5 J/g and the melting point, T f , is 19 • C [49]. The thermophysical properties of the constituents are shown in Table 2. In our modeling approach, we considered 2D unit cells of the material including the PCM and the cementitious matrix, which were explicitly modeled at different volume fractions. For the FEM discretization, we considered isoparametric quadrilateral elements. Figure 2 shows three different meshes. The material properties considered for each phase were those from Table 2.
Heat transfer is described as an elliptical partial differential equation that reads: where ρ is the density, c p the specific heat, T is the temperature field, t is time, k is the thermal conductivity, and q is the volumetric heat source.
In this work, we account for the non-linear effect of the phase transition following the effective specific heat approach [51,52]: where c p,s and c p,l are the specific heat in the solid and liquid phases, respectively, h f is the enthalphy of fusion, T f is the phase-change temperature, and ∆T f the temperature window. In order to solve an initial boundary value problem (IBVP), apart from the partial differential equations involved, we need to define the solution domain (i.e., the unit cell and time) and initial and boundary conditions. In our simulations we applied Dirichlet boundary conditions on the left side of the unit cell and evaluated the temperature evolution at different locations ( Figure 3).
The heat transfer problem can be discretized using the FEM [53], yielding the following matrix form: where C th and K th are the heat capacity and thermal conductivity matrices, respectively, T is the temperature vector, and q is the flux vector. This is a time-dependent problem; therefore, we can discretize in time using Equation (3) with a time increment ∆t, such that at an instant of time t j , the system reads: We make use of the Crank-Nicolson time integration scheme [53] to determine the temperature at an instant of time t j+1 : Since this problem is non-linear, we must update the capacity matrix at each time step, i.e., C th = C th (T j ), which gathers the phase-change phenomenon, as described above.

Upscaling Procedure
The unit cell model is heterogeneous, since it accounts for two different phases. However, for the component design, we require material properties of the mix. For this purpose, we use an upscaling homogenization approach [54] for determining the macroscopic thermophysical properties, namely, k and ρc p .
The effective conductivity, k e f f , of a two-phase system with spherical inclusions can be obtained as [55]: where k m and k p are the thermal conductivity of the matrix and particle, respectively, and φ is the volume fraction of particles. In our case, V f = φ, so we keep using φ for the mathematical description of the problem in the following. The effective volumetric heat capacity, ρc e f f p , is evaluated using mixing rule [52]: where ρ m and ρ p are the densities of the matrix and particle, respectively. Similarly, c m p and c p p are the specific heats of the matrix and particle, respectively.
Using the definitions above, it is possible to set the homogenized phase-change (i.e., temperature-dependent) model as: where ρc e f f p,s and ρc e f f p,l are the effective volumetric heat capacity of the solid and liquid phases, respectively, ρ p is the density of the particle, and h f and T f are the heat of fusion and melting point of the PCM.

Component Design
The new building envelope proposed herein is made of a PCM-based cementitious matrix, described previously, with air cells, yielding a honeycomb-like structure. In the first place, tailor-made molds were fabricated using 3D printing technology. Then, different brick configurations were cast and tested experimentally using a hot plate setup. Finally, the designs were validated using a macroscopic FEM model. Moreover, the simulation framework was also validated experimentally so that it can be used in the analysis of different designs.

Experimental Methods
The building envelopes follow a honeycomb structure so that large amounts of air could be hosted within the blocks while preserving a compact and robust structure. Due to the complex geometry of the blocks, new tailor-made molds were fabricated using 3D printing technology [56] (Ultimaker 2 + 3D) with acylonitrile butadiene styrene (ABS) and polylactic acid (PLA). Namely, we used ABS for the large air cells design, and we used PLA for the other two. The advantage of PLA over ABS is that the former is biodegradable and allows a higher print speed [57]. In any case, the choice of the material does not affect the results of this research. Three blocks with size 100 × 52 × 20 mm 3 and different air cell sizes were designed using CAD techniques ( Figure 4). Once the bricks were printed, these were filled with a mix of two components of PS8510 silicon rubber, which takes two hours to harden. In such a way, molds of the printed bricks were realized in a flexible material, which is better for the demolding process ( Figure 5a). Another option is to directly print the molds in 3D, but this may result in lower precision [58]. Finally, three types of honeycomb blocks (i.e., small, medium, and large cells), as shown in Figure 5b, were cast using the material mix described previously. Ordinary cement paste (CEM I 42.5, w/c = 0.45) included encapsulated organic paraffin wax with different volume fractions (V f = 0%, 10%, and 30%). Table 3 shows the different configurations considered in our analysis. To analyze the thermal behavior of the new blocks, a hot plate setup ( Figure 6) was used along with a FLIR A320 thermal camera and ThermaCAM Researcher Pro software. All the specimens were insulated on the free boundaries (i.e., back, top, and sides) except for the front, which was left open so as to evaluate the temperature evolution with the thermal camera. The insulation was made of styrofoam, sponges, and duct tape. The hot plate was set to 100 • C, and the room temperature was 22.5 ± 1 • C. The tests had a total duration of 3.5 h.
The evolution of the temperature was obtained for three sample heights, namely, 30%, 60%, and 90%, to see the spatial gradient. Temperature was evaluated every 9 min until the steady-state conditions were accomplished.

Brick FEM Model
At the component level, we used the homogenized phase-change FEM formulation presented in Section 2.1.2 for the non-linear unsteady heat transfer problem. For the FEM model, we used isoparametric quadrilateral elements, and the boundary conditions reflected those of the hot plate setup.
As shown in Figure 7, the contact of the bottom side of the brick with the hot plate was not perfect, and internal convective losses due to the surface roughness appeared. In order to model this, we used thermal contact resistance elements (i.e., interfacial conductance) that were calibrated with part of the experimental results. The temperature of the source was set to that of the hot plate, 100 • C. Air trapped in the cells would also generate convective losses that would affect the heat transfer process in the blocks. These were taken into account in our model by using film conditions: where q bottom,j is the imposed flux at the bottom nodes and instant t j , h c is the interfacial conductance, A bottom is the contact surface, T plate is the temperature of the hot plate, T bottom,j is the temperature of the bottom surface of the brick, q cells,j is the imposed flux on the inner surface of the cells, A cells is the inner surface of the cells, T air is the air temperature at the cells, and T cells,j is the temperature at the inner surface of the cells.
The time evolution of the temperature distribution was analyzed by placing three probes in the domain at the center of the horizontal axis, and three different heights: 30%, 60%, and 90%-those recorded in the experimental setup.
We extended our analysis of the component so as to evaluate the mechanical stresses promoted by the temperature gradients. Thus, according to the constitutive equation in solid mechanics: where σ is the stress tensor, D the consitutive tensor, and ε the elastic strain tensor. The elastic strain tensor can be evaluated as: where ε is the total strain obtained from the displacement field gradient, and ε th is the thermal strain tensor, which in isotropic solids is: where α is the thermal expansion coefficient; ∆T is the temperature jump defined as the difference between the temperature field obtained from the heat transfer analysis, T, and the reference temperature, T r , i.e., ∆T = T − T r ; and I is the identity tensor. Regarding the FEM model, we used isoparametric quadrilateral elements for the space discretization, and the system was solved at every instant of time t j : where K is the stiffness matrix, u j the displacement vector, and f j the force vector, which is defined as: with f m,j being the mechanical force vector and f th,j the thermal force vector.
For an element el and isotropic expansion, the element thermal force vector, f el th,j , reads: B is the matrix containing first-derivatives of the shape functions, and D is the constitutive matrix (in our case, plane stress problem). Finally, we solved the mechanical problem for each t j using the preconditioned conjugate gradient method.

Material Thermal Characterization
In the first place, we performed heat transfer transient simulations at the unit cell scale to characterize the thermal behavior of the two materials designs. For this purpose, we defined three unit cells with the following volume fractions of PCM: V f = 0%, 10%, and 30%. The size of the unit cell was 10 µm, and the capsule sizes were 3.5 and 6 µm for the 10 and 30% contents, respectively, which are in the ranges of similar compositions [59][60][61]. The FEM mesh consisted of approximately 2000 elements.
The unit cells were initially set to T(x, 0) = 15 • C, and the temperature at the left edge, x = 0, identified with the label "a" in Figure 3, was suddenly set to T(0, y, t) = 30 • C. The non-linear heat transfer problem was solved for the first 200 s. As observed in Figure 8, at the microscale it took approximately 36 s to complete the phase transition, which in our paraffin wax takes place at T f = 19 • C. We have considered a temperature window of ∆T f = 2 • C and the thermophysical properties in Table 2.  Figure 9 shows the evolution of the temperature of the points labeled as a and e as per Figure 3, i.e., the left and right edges of the unit cell. We used two types of model: (i) a heterogeneous model accounting explicitly for the phase-change phenomenon in the PCM particles (results are shown as continuous lines in Figure 9), and (ii) a homogeneous model using the homogenization approach presented in the previous section (results are shown as dotted lines in Figure 9). We can observe how the temperature rises during the first 9 s from 15 • C to 18 • C, where the phase transition starts taking place in our approach (we use a temperature window of ±1 • C with respect to the melting point (T f = 19 • C). Thus, we observed a change in the slope of the curves when the phase transition, solid to liquid, started.
In the case of no PCMs (V f = 0%), it can be seen that there was no change in slope, and the temperature at the right edge reached that of the one imposed at the left edge at 200 s. This means that steady-state conditions were reached. As we increased V f , we observed that the material required more time to achieve these conditions. Thus, the addition of PCM delays the rising of temperature, and we also confirmed that the higher the PCM amount, the higher the thermal inertia exhibited by the composite material. In fact, with the heterogeneous model, we measured 1.8× delay until reaching 95% of the steady-state conditions for V f = 10%, and this number increased up to 3.4× in the case of V f = 30%. We also observed some differences between the homogeneous model, in which one phase with homogenized thermophysical properties was modeled, and the heterogeneous one. These differences may have come from the fact that the homogenized expressions correspond to a 3D case (inclusions were treated as spheres), whereas the heterogeneous model considered inclusions as disks. In any case, the same meshes were used for both models. With the the homogenized model being more conservative in terms of thermal management capabilities, we used it in the brick simulation.

Brick Characterization
Once the new mix (i.e., PCM and cement) had been characterized, we proceeded with the characterization of the building envelopes. For this reason, we tested different designs (summarized in Table 3) using the hot plate setup. We also performed FEM simulations of the bricks to better understand the effects of V f and φ air on the global insulating properties of the bricks. As pointed out in Section 2.2.1, thermographic analysis was used to measure the temperature evolution within the brick. Figure 10 shows the temperature distribution obtained through thermographic analysis and the corresponding estimates using the FEM model with homogenized properties. As a first step in our modeling tasks, we needed to calibrate the thermal conductance of the brick-plate interface. In Figure 11 we show the experimental results (in dotted lines) and the equivalent simulation results (in continuous lines) of the large air cell model, which was the one used for the calibration, since we have it for the three different grades of PCM. From the results, we can observe that a conductance of h c = 200 W/m 2 K provides a good fit. This value of h c is in agreement with that observed for concrete-steel interfaces [62]. We considered an air film coefficient of h air = 10 W/m 2 K.   Figure 11. Numerical vs. experimental temperature evolution for V f = 0% and large air cells at different heights of the brick.

Effect of V f
In Figure 12 we show the temperature evolution during ca. 2.5 h, for the large air cell specimens with V f = 10% (Figure 12a) and V f = 30% (Figure 12b). In the first place, we can observe that the increase in PCM content provided higher thermal inertia and insulating properties. In fact, during the first 3000 s, the temperature rise in the mix with higher V f increased at a lower pace than that of V f = 10%. The steady-state conditions were reached after close to 4000 s.
If we look at the 90% height, which is the most interesting line, for it is closest to the indoor one, we can observe again that after 1000 s, the difference in temperature was about 2.5 • C between both mixes. The steady-state conditions were also milder in the case of V f = 30%.

Effect of Air Cells
We now study the effect of air cells in the insulating properties of the bricks. In Figure 13a it can be seen that at 60% of the height of the samples without PCM, the temperature difference after 6000 s was greatest between the samples with the small air cells and the large air cells. The sample with the large air cells heated up the slowest, as expected, since the volume of the air cells covered 44% of the total volume of the cementitious sample, whereas the small air cells and the medium air cells only covered 30% and 24% of the total volume, respectively. The fact that the sample with the mediumsized air cells heated up slower than the one with the small air cells does not seem to be in line with the literature. The same seems to happen in Figure 10, where the average temperature at 90% of the height of the samples is shown. The sample with small air cells heated up faster than the sample with medium air cells, which is not in proportion with the volume fraction of air within the different samples. Just like in the case of the PCM-containing samples, the graphs for 60% and 90% height are very similar in terms of shape and temperature difference between the samples.
In Figure 13b, it can be observed that larger hole specimens (L) presented a temperature decrease of at least 5 • C with respect to small hole specimens (S) when no PCM was added, and up to 7.5 • C with V f = 30%. Large, V f =10% Large, V f =30% Figure 13. Experimental temperature evolution of: (a) V f = 0% and (b) All specimens at 90%.

Thermomechanical Assessment
Finally, we used the temperature field obtained from the heat transfer simulations to run a sequential FEM of the mechanical problem. We considered an elastic modulus E = 20 GPa, as observed by [49], a Poisson's ratio ν = 0.2, and thermal expansion coefficient α = 10 −5 1/K. We ran the simulation for the first 3000 s, where the most important temperature gradients appeared.
In Figure 14, we show the total deformation of the brick, for large air cells and V f = 10%. One can see the free expansion deformation mode with the highest values at the bottom tips. Figure 15 shows the maximum (Figure 15a) and minimum (Figure 15b) principal stresses at the same instant of time. As expected from the deformation field, the highest maximum principal stress appeared at tips of the bottom surface. It must be remarked that we considered sliding boundary conditions for the bottom surface. Apart from the stress concentration due to boundary conditions, we observed well-developed stress bands below the tensile strength (ca. 2 MPa). Regarding the minimum principal stresses, these are well below the compressive strength as measured in [49]. This is in agreement with the experimental observations, where no cracks were identified.

Conclusions
The design of new insulating envelopes is a direct route towards energy efficient buildings. The use of advanced numerical and manufacturing techniques may lead to important advances in this sense. In this work we combined in a seamless fashion new materials (e.g., PCMs), manufacturing (e.g., 3D printing), and simulation techniques (e.g., finite element method) in the design of new insulating components with enhanced properties. Namely, we presented the so-called HEXCEM component, which has important benefits, such as improved thermal insulation, reduced energy consumption and GHG emissions, reduced specific weight, and applicability (it could be also used in new buildings as well as in the rehabilitation of older buildings).
For this purpose, a new synthesis approach based on additive manufacturing of tailormade molds and microencapsulated PCM has been successfully implemented. As we showed in this work, complex geometries can be easily accomplished.
In order to reinforce the design process for both the material mix and the component, a two-scale (material and component) simulation strategy has been successfully developed. For that, the numerical simulations carried out with our methodology showed excellent agreement with the experimental tests. Moreover, it can be used to explore new designs for different applications before actually synthesizing them.
From the analyses, we found that microencapsulated PCM in cement, within the phase-transition temperature, may lead to 1.8× delay to 95% steady-state conditions in the case of V f = 10%, and 3.4× for V f = 30%. On the other hand, larger hole specimens (L) present a temperature decrease of at least 5 • C with respect to small hole specimens (S) when no PCM is added, and up to 7.5 • C with V f = 30%.
With respect to the thermomechanical behavior, the simulations of the bricks showed no thermal cracking. This is also in agreement with the experimental observations. This is important for ensuring mechanical integrity.