{'fraction_non_alphanumeric': 0.03614288393024564, 'fraction_numerical': 0.01898556159759985, 'mean_word_length': 5.160265665607854, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Microstructure-guided numerical simulation to evaluate the influence of phase change materials (PCMs) on the freeze-thaw response of concrete pavements

: The use of phase change materials in infrastructure has gained significant attention in the recent years owing to their robust thermal performance. This study implements a numerical simulation framework using finite element analysis to evaluate the influence of Phase Change Materials (PCMs) on the thermal response of concrete pavements in geographical regions with significant winter weather conditions. The analysis is carried out at different length scales. The latent-heat associated with different PCMs is efficiently incorporated into the simulation framework. Besides, the numerical simulation framework employs continuum damage mechanics to evaluate the influence of PCMs on the freeze-thaw induced damage in concretes. The simulations show significant reductions in the freeze-thaw induced damage when PCMs are incorporated in concrete. The numerical simulation framework, developed here, provides efficient means of optimizing the material design of such durable PCM-incorporated concretes for pavements by tailoring the composition and material microstructure to maximize performance.


INTRODUCTION
Use of Phase change materials (PCMs) in cementitious composites for enhanced energy efficiency is an active research area [1][2][3][4]. PCMs have been shown to be advantageous in structural concretes due to their ability to store and release heat [2,[4][5][6]. A large number of studies have focused on the use of PCMs towards indoor thermal comfort in buildings and enhancement of building energy efficiency [6][7][8][9]. Besides, PCMs have been shown to reduce early age hydration heat release and the early-age thermal cracking in large concrete structures and pavements [10][11][12]. In addition to the above-mentioned applications, availability of these PCMs with a wide range of transition temperatures [1,13] opens up various innovative avenues towards beneficial utilization of these materials in infrastructure materials. For example, application of PCMs towards reduction in freeze-thaw cycles in concretes pavements has been proposed earlier [4,13].
Deicing salts are commonly used to remove snow and ice from the pavement surfaces [4,[13][14][15][16][17]. But, such process is labor-intensive and causes damage to the pavements [4,13,18,19]. Alternatively, incorporation of PCMs in concretes has been suggested [1,4,20] to solve this issue by the means of high latent heat of fusion of PCMs. A recent experimental study [4] has shown significant reduction in snow formation when PCMs are incorporated in concrete. Another recent study confirmed improved frost resistance of the concrete containing PCMs [20]. Numerical studies involving freeze-thaw performance of cementitious materials are also implemented using one dimensional finite difference technique [10,13,20]. This paper focuses on microstructure-guided numerical simulation using finite element analysis to evaluate the influence of microencapsulated paraffinic PCMs with transition temperature around 5℃ on the freezethaw response of macro-scale concrete pavements. The numerical simulation framework, presented in this study, performs numerical homogenization using finite element analysis at multiple length scales. The effective properties, obtained using multi-scale numerical homogenization, are employed in a macro-scale FE model of concrete pavement to obtain pavement temperatures under the effects of the ambient temperature and solar radiation whilst considering conduction, convection and radiation modes of heat transfer with the subgrade and the atmosphere. The pavement temperatures, thus obtained, are used to quantify frost damage at multiple length scales using continuum damage mechanics [21][22][23][24] that can capture the progressive freeze-thaw induced damage. Thus, this study integrates the material microstructure and thermo-mechanical properties of the constituent phases to predict the freeze-thaw response of macroscale concrete pavement. The numerical simulation framework, presented in this study, provides efficient design strategies to tailor the composition and microstructure towards development of concretes with improved frost-damage performance without compromising its strength. Hence, the methodology developed herein would be extremely useful for pavements in areas with severe winter weather conditions.

NUMERICAL SIMULATION TO ELUCIDATE THE INFLUENCE OF PCM ON THERMAL RESPONSE OF CONCRETE PAVEMENTS
In this paper, a thermal analysis is carried out at multiple length scales with a view to evaluate the influence of PCMs in concrete pavements. Here two different material mixtures are considered. While the control concrete contains 70% aggregates (40% coarse aggregate and 30% fine aggregate) by volume, the PCMincorporated concrete contains 64% aggregates (40% coarse aggregate and 24% fine aggregate) and 6% microencapsulated paraffinic PCMs (as 20% volumetric sand-replacement) by volume. Such dosage of PCM has been shown to provide equivalent mechanical performance as compared to conventional mortars [25]. In this study the dosage of PCMs is fixed at 20% by volume sand substitute since incorporation of higher volume fraction of PCMs has been shown to result in decrease in compressive and flexural strength in a previous study [25]. PCM considered here has phase transition temperature of 5.1℃. It consists of a paraffin core encapsulated by melamine formaldehyde shell as detailed in [25].
The forthcoming sub-sections describe numerical homogenization for effective property-prediction and application of effective properties towards thermal response evaluation of a macro-scale pavement section with a view to highlight the influence of PCMs in concrete.

Multiscale Numerical Homogenization for Prediction of Effective Thermal Properties:
Here, a microstructure guided numerical simulation framework is employed at multiple length scales to compute the effective properties of both the concretes. At every length scale, a representative volume element (RVE) is generated with the inclusion phases. Intrinsic thermal and mechanical properties are then assigned to the component phases in the RVE. Periodic boundary conditions [26,27] are applied. The RVE is meshed and FE analyses are carried out. A post processing module generates volume averaged responses that are computed to yield effective properties. The basic framework for the procedure is illustrated in The following sub-sections present the multiscale numerical homogenization procedure to obtain the effective properties of both the concretes.

Generation of Microstructure and Boundary Conditions:
3D representative volume elements (RVEs) are generated using the Lubachhevsky-Stillinger algorithm [28][29][30][31]. A hard particle contact model is employed, while ensuring there is no particle overlap. The formulations regarding the generation of the RVEs are adequately described in [26,27,32]. The meshed microstructure is implemented in a Python script and the analysis is performed in ABAQUS™ solver.
Periodic boundary conditions (PBC) [33][34][35][36] are applied on the RVE. Periodicity provides higher computational efficiency in smaller analysis domain and thus facilitates faster convergence [27]. Periodic boundary conditions (PBC) ensure continuity of temperature and heat flux at the boarder of neighboring unit cells. PBCs have been applied successfully towards prediction of effective properties of heterogeneous composites [26] and are detailed adequately in [33][34][35][36]. Figure 2 presents the microstructures generated for multi-scale homogenization. Step I: air voids in HCP; (b-2) PCM in homogenized matrix from step-I ;(b-3) sand in homogenized matrix from step-II ; (b-4) coarse aggregate in homogenized matrix from step-III A three step homogenization procedure is followed for the case of control concrete. The first step homogenizes air voids with hardened cement paste (Figure 2 (a-1)) whereas the second step homogenizes sand into the homogenized matrix obtained from first step (Figure 2 (a-2)). The final step incorporates coarse aggregates into the homogenized matrix obtained from second step (Figure 2 (a-3)) and thus, effective property of concrete is computed.
While a three-step homogenization is employed for control concrete, an additional step involving homogenization of PCMs is adopted for PCM-modified concrete. The first step of homogenization remains the same as control concrete case. In the second step, microencapsulated PCMs are incorporated into the homogenized matrix obtained from the first step as shown in Figure 2 (b-2). The third and fourth steps homogenizes sand ( Figure 2 (b-3)) and coarse aggregate ( Figure 2 (b-4)), respectively, into the homogenized matrix obtained from the previous step. These homogenizations are performed here at multiple steps in order to reduce the computational cost while maintaining computational efficiency [37].
The relative edge lengths of RVEs with respect to inclusion sizes are computed (Shown in Figure 2) based on sensitivity study. Similar relative sizes of RVEs with respect to inclusion sizes were obtained elsewhere [26,27]. The size distribution of inclusions (air voids, PCMs, sand and coarse aggregates) are adopted from [22,25,26]. The absolute volume fractions of component phases are shown in

Effective Property Computation
Once the RVE is generated using the Python script, it is meshed and subjected to FE analysis under externally applied unit temperature gradient along the X-direction. The FE analyses are carried out in ABAQUS™ solver. A MATLAB© script operates on the result files containing the volume of the elements and the responses to obtain the effective properties as explained before. The intrinsic thermal properties of the constituent phases of both the concretes are shown in the Table 2. Table 2. Thermal conductivity ( ) [26], specific heat capacity (C) [5] and latent heat L [4,13]   Step-IV: coarse aggregate in homogenized matrix from step-III: Homogenized : 1.78 W/m K. All the heat flux values correspond to externally applied unit temperature gradient (1℃/ ) in X-direction.
A post processing module (MATLAB© subroutine) computes the volume average of the heat flux. The volume-averaged heat flux is used to obtain the effective thermal conductivity using Equation 1 [26].
where is the effective thermal conductivity, ̅ is the homogenized (volume averaged of all elements) heat flux in RVE, L is the edge length of RVE. The temperatures on the two opposite faces of the unit cell along the X-direction are denoted by =0 and = . Thus, difference. The aforementioned 3D numerical homogenization technique for effective thermal conductivities has been successfully validated in a previous study [26] against experimental thermal conductivities obtained using a guarded hot plate apparatus in accordance with ASTM C177-13 [38]. The effective thermal conductivities of a multi-inclusion porous systems is often predicted on the basis of a microstructure guided study [26,39] which have shown substantial agreement with experimental results.
In order to obtain the effective volumetric heat capacity, the reaction heat flux is obtained at the nodes under applied unit temperature gradient along X-direction. The ratio of the volume-averaged reaction heat flux to the temperature difference yields the effective volumetric heat capacity ( ) , [40]. For the case of PCM incorporated concrete, a similar procedure is followed to obtain ( ) , that has PCM in the solid state. However, at the transition temperature of the PCMs, the latent heat released by the PCM during phase transformation is appended to the volumetric heat capacity as shown in Equation 2 [5].
where ( ) , is the effective volumetric heat capacity at transition temperature, ( ) , is the effective volumetric heat capacity for solid phase PCM, is volume fraction of PCM, is density of PCM, ℎ is latent heat of fusion of PCMs and ∆ is the phase change temperature window. It is to be noted that all the calculations involved herein consider PCMs in solid state only, for the sake of simplicity.
The effective properties obtained using the above-mentioned framework for both control and PCMmodified concretes are presented in Table 3. The values obtained for control concrete are in line with experimental values reported in literature [5,13]. For the cases of PCM inclusions, the values match closely with the ones reported in [5,26]. The effective thermal properties, thus obtained, are used in macro-scale simulation of pavement in the forthcoming section to elucidate the influence of microstructural modifications imparted by PCMs on the macro-scale structures.

Model Geometry and Boundary Conditions:
A three-dimensional pavement section is modeled in FE commercial software ABAQUS™ as shown in  define the heat transfer interactions. Convection is defined by the product of surface convection coefficient (ℎ 0 ) and the difference of the pavement surface temperature ( ( , )) and ambient temperature ( 0 ) [44].
The value of ℎ 0 is a function of wind velocity which can be assumed to be 26 W/m 2 K for the air-concrete interface at low wind speeds [45] and is applied as a surface film coefficient in the model that  Table 3), conductivity of subgrade is adopted as 1.75 W/(m.K) [47]. Thus, as the pavement in the model cools down in absence of solar radiation and plunging ambient temperature, the PCMs having a volume fraction of 6% (Vf), release 150 kJ/kg (hsf) at 5.1℃ (Tpc). Density of 300 kg/m 3 for PCMs is adopted here [25]. This study incorporates a phase transformation temperature window (∆ ) of 1℃ which is in line with values reported elsewhere [5]. At the transition temperature of the PCMs, the latent heat released by the PCM during phase transformation is added to the volumetric heat capacity (Equation 2), obtained using numerical homogenization as explained before.

Finite Element Simulation of Macro-scale Response:
A MATLAB© script extracts the data from the weather data files and feeds it into the ABAQUS™ model.  The ambient temperature is also depicted in Figure 5 for reference. As can be seen from the figure, the pavement surface temperature for both concrete systems is higher than the peak ambient temperature which can be attributed to the solar radiation that heats up the surface. Similar observations were noted elsewhere [10,13]. Figure 5 also clearly shows that temperature peak of control concrete is higher than the temperature peak of the PCM-incorporated concrete when the temperature is greater than the phase transition temperature of 5.1℃. The reduced temperature peak obtained in the PCM incorporated concrete substantiates the effective solar energy storage capacity of PCMs which would otherwise be wasted in raising the surface temperature as obtained in the control concrete. The effectiveness of the PCM incorporated concrete lies in its ability to store the energy during its phase change as a latent heat in the  [13]. In order to shed more light, Figure 6 shows the pavement temperature profiles for January 7, 2016.  This is attributed to the heat released by PCMs during phase-transformation thus helping to keep the pavement warmer. Considering a water-freezing temperature of -6.1℃ in concrete [13,48], the control concrete pavement is expected to freeze at this minimum ambient temperature of on January 7, 2016 whereas the PCM-incorporated concrete does not freeze. Overall, considering a freezing and thawing temperature of -6.1℃ and 0℃ respectively in concrete [13,48], the control concrete experiences 25 freezethaw cycles, whereas the PCM-incorporated concrete experiences 13 number of freeze-thaw cycles during December 2015 -February 2016 (Surface temperatures reported in Figure 5) thereby reinforcing the idea that PCMs have the potential to enhance the service life of pavements. In order to shed more light into the effectiveness of the PCMs, simulation of frost-damage in both the concretes is elucidated in the forthcoming section. during the freezing cycle leads to a volume increase of 9% in the micro-pores [49]. The damage in the cementitious matrix of the REA is implemented in ABAQUS™ with the help of a user defined subroutine UMAT that implements the damage laws. A post processing module computes the degraded Young's modulus from the results obtained by ABAQUS™ solver. The analysis procedure and the results are presented as follows:

Microstructures and Material Properties:
The damage analysis is carried out at the mortar scale. Here, 2D REAs are chosen to reduce computational demand as a trade-off between computational efficiency and demand [50]. The REAs are generated using Lubachevsky-Stillinger algorithm [26,[28][29][30][31]. The edge lengths of REAs are chosen as five times the mean diameter of inclusions based on a sensitivity study. Similar RVE sizes relative to inclusion sizes are implemented successfully in [27,32]. Figure 7(a) shows the control mortar microstructure containing 40% hardened cement paste with 10% voids and 50% sand by volume. The REA for PCM-incorporated mortar is shown in Figure 7(b) which contains 20% PCM as sand-replacement. The PCM-modified mortar microstructure contains 40% hardened cement paste with 10% voids, 40% sand and 10% PCM (20% sand replacement) by volume. The size distribution of PCM and sand are adopted from [25]. The analysis requires the material properties for the matrix and inclusion phases which are shown in Table   4. The poison's ratio for all components except PCMs is taken as 0.2. A poison's ratio of 0.4 is adopted for PCMs [25]. Table 4. Input material properties for numerical simulation of damage [22,25,37]

Numerical Simulation of Isotropic Damage in Concretes:
The damage analysis is implemented in python language program which is later used in ABAQUS TM solver to obtain results. A user-defined subroutine in ABAQUS TM applies the volume expansion (due to freezing) in the voids embedded in the cement paste matrix. The resulting mechanical response is captured by the framework that applies the damage law. The mechanical properties are used from Table 4 while the size distribution is similar to what has been implemented earlier in the study. Finally, a post processing module computes the stress-strain response after every freeze-thaw cycle and yields the degraded material properties of the mortar which is scaled up to the concrete level and validated with experimental results reported in [49].
Periodic boundary conditions [33][34][35][36] are implemented in the REAs as described earlier. The framework assumes fully saturated voids. The analysis framework imports the surface temperatures of the pavement from the macro scale analysis (for the month of February 2016 as shown in Figure 5(c)) and applies volumetric expansion in the voids as soon as the freezing temperature of -6.1℃ [13,51] is reached. The volumetric expansion is implemented using a user-defined subroutine in ABAQUS TM that incorporates volumetric thermal expansion of voids by means of artificial coefficient of thermal expansion of the void so as to obtain 9% volumetric expansion. The volume expansion is released when the average nodal temperature increases above 0℃ (Thawing temperature [13,52]). Such freeze-thaw cycle induce damage around the voids which propagates with increasing number of freeze-thaw cycles.
The progressive damage in HCP matrix due to freeze-thaw cycles is implemented by user defined subroutine UMAT [53][54][55] in ABAQUS™. The progressive degradation of both the mortars due to frostdamage is modelled here by continuum damage mechanics while considering the heterogeneity in these materials [21][22][23]. Assuming isotropic stiffness degradation, the damage variable, D is given as [21][22][23][24]56]: Where the effective stress tensor is denoted by , ∁ denotes fourth order tensor of elasticity and is the strain tensor. The value of damage, D ranges from 0 (undamaged) to 1 (completely damaged). The damage rate denoted by ̇ assumes only zero or positive values which can be explained by its proportionality with the damage energy release rate [22].
In order to characterize the damage in the cement paste matrix at the scale of the mortar (Figures 7(a) and , an isotropic damage model is used that provides a damage evolution law as shown in Equation 5 [57].
̇= ( ( ) ⁄ ) [5] Where is the Lagrange multiplier which when multiplied with the partial derivative of the constraint ( ) (as mentioned in Equation 6) with respect to the (elastic energy) gives the rate of scalar damage ̇. The constraint ( ) is defined as follows: ( ) is constrained to be ≤ 0. The damage evolution law (Equation 5) is discretized by an implicit Euler scheme and it yields the instantaneous damage parameter . The fitting parameters a, b and c (in Equation   6) for cement paste, adopted from [57], are 0.488, 0.1 and 3.0 respectively. For the PCM incorporated mortar, same damage parameter values are implemented since the matrix is HCP for both the cases.   presented in this paper, provides an efficient means of designing concretes by tweaking the material at different length scales to maximize performance.